complex_matrix_market_demo/

complex_matrix_market_demo.rs

//! Complex Matrix Market benchmark/demo for GMRES/FGMRES on generated and file-backed
//! complex systems. Run it with MPI enabled when comparing replicated correctness
//! checks against distributed scalability behavior.
//!
//! Correctness sweep (replicated operator/preconditioner checks, long iteration
//! budget, and stricter modified Gram-Schmidt reorthogonalization):
//! cargo mpirun -n 4 --example complex_matrix_market_demo --features=complex,mpi,mpi_examples -- --mode correctness --dist-policy off --maxits 2000 --restarts 150,250,324 --pcs replicated-full-ilu0,none,jacobi,block-jacobi-ilu0 --fgmres-orthog modified --fgmres-reorth always
//!
//! Scalability sweep (automatic distributed policy with lightweight local-block
//! preconditioners):
//! cargo mpirun -n 4 --example complex_matrix_market_demo --features=complex,mpi,mpi_examples -- --mode scalability --dist-policy auto --maxits 1000 --restarts 50,100 --pcs none,jacobi,block-jacobi-ilu0
//!
//! Optional MPI correctness sweep with overlap ILU rows. `--include-overlap-pcs`
//! keeps the zero-overlap block-Jacobi ILU baseline, then adds
//! `asm-ilu0-overlap1`, `ras-ilu0-overlap1`, and `ras-iluk1-overlap1`;
//! correctness rows include `GlobalOK`/`XOK` verdicts before users compare speed.
//! cargo mpirun -n 4 --example complex_matrix_market_demo --features=complex,mpi,mpi_examples -- --mode correctness --dist-policy auto --include-overlap-pcs --warmup-runs 0 --measured-runs 1

#![cfg_attr(not(feature = "complex"), allow(dead_code))]

#[cfg(not(feature = "complex"))]
fn main() {
    eprintln!(
        "This example requires the `complex` feature. \\nre-run with `cargo run --features complex --example complex_matrix_market_demo`."
    );
}

#[cfg(feature = "complex")]
use kryst::error::KError;

#[cfg(feature = "complex")]
mod complex_demo {
    use std::env;
    use std::fs::File;
    use std::io::Write;
    use std::path::{Path, PathBuf};
    use std::sync::Arc;
    use std::time::Instant;

    use super::KError;
    use kryst::algebra::bridge::BridgeScratch;
    use kryst::algebra::prelude::*;
    use kryst::context::ksp_context::{ReorthPolicy, Workspace};
    use kryst::matrix::DistCsrOp;
    use kryst::matrix::dist_csr::DistributedPlanDiagnostics;
    use kryst::matrix::sparse::CsrMatrix as SparseCsrMatrix;
    use kryst::ops::klinop::KLinOp;
    use kryst::ops::kpc::KPreconditioner;
    use kryst::parallel::{Comm, UniverseComm};
    use kryst::preconditioner::PcSide;
    use kryst::preconditioner::Preconditioner;
    use kryst::preconditioner::ilu_csr::{
        IluCsr, IluCsrConfig, IluKind, ReorderingKind, ReorderingOptions,
    };
    use kryst::preconditioner::jacobi::Jacobi;
    use kryst::preconditioner::overlap_ilu::{OverlapIluPc, OverlapRestriction};
    use kryst::solver::fgmres::{
        FgmresSolver, FgmresStagnationPolicy, FgmresVariant, OrthogMethod, ResidualCheckPolicy,
    };
    use kryst::solver::gmres::{GmresOrthog, GmresSolver, GmresVariant as PlainGmresVariant};
    use kryst::solver::{LinearSolver, MonitorAction, MonitorCallback};
    use kryst::utils::convergence::ConvergedReason;
    use kryst::utils::matrix_market::read_matrix_market;

    #[cfg(feature = "mpi")]
    use kryst::parallel::MpiComm;
    #[cfg(not(feature = "mpi"))]
    use kryst::parallel::NoComm;

    pub fn run() -> Result<(), KError> {
        #[cfg(feature = "logging")]
        let _ = env_logger::try_init();

        #[cfg(feature = "mpi")]
        let comm = UniverseComm::Mpi(Arc::new(MpiComm::new()));
        #[cfg(not(feature = "mpi"))]
        let comm = UniverseComm::NoComm(NoComm);

        let rank = comm.rank();
        let size = comm.size();
        let config = BenchmarkConfig::from_env_args()?;
        #[cfg(feature = "mpi")]
        let is_parallel = matches!(comm, UniverseComm::Mpi(_)) && size > 1;
        #[cfg(not(feature = "mpi"))]
        let is_parallel = false;

        if rank == 0 {
            println!("Complex Matrix Market demo (GMRES/FGMRES)");
            println!("Run mode: {}", config.run_mode.label());
            println!(
                "Parallel backend: {}",
                if is_parallel {
                    format!("MPI ({} ranks)", size)
                } else {
                    "serial".into()
                }
            );
            println!(
                "Benchmark runs: warmup={}, measured={}",
                config.warmup_runs, config.measured_runs
            );
            println!(
                "Mode profile: {}",
                match config.run_mode {
                    RunMode::Correctness => {
                        "small matrices, explicit true residuals, optional replicated operator checks"
                    }
                    RunMode::Scalability => {
                        "distributed operator + local-block PC only, global norms only, lightweight reporting"
                    }
                    RunMode::Robustness => {
                        "robustness stress: fallback-enabled stagnation handling and restart-heavy behavior"
                    }
                }
            );
            println!(
                "Stagnation fallback: {} (min_inner_before_fallback={})",
                if config.run_mode == RunMode::Correctness || !config.allow_stagnation_fallback {
                    "disabled"
                } else {
                    "enabled"
                },
                config.min_inner_before_fallback
            );
            println!("KSP default: {}", config.ksp.label());
            println!("GMRES/FGMRES haptol: {:.3e}", config.fgmres_haptol);
            println!(
                "Complex ILU reordering: kind={:?}, symmetric={}",
                config.ilu_reordering.kind, config.ilu_reordering.symmetric
            );
            if config.run_mode == RunMode::Correctness {
                println!(
                    "Replicated check marker: {}",
                    if config.mark_replicated_check {
                        "enabled (metadata-only)"
                    } else {
                        "disabled"
                    }
                );
                println!(
                    "Overlap PC preset: {}",
                    if should_include_overlap_pcs(&config, size) {
                        "enabled (MPI correctness only; GlobalOK/XOK verdict columns shown)"
                    } else if config.include_overlap_pcs {
                        "requested but inactive (requires MPI correctness with >1 rank)"
                    } else {
                        "disabled"
                    }
                );
            }
            if config.residual_history {
                println!(
                    "Residual history: enabled{}",
                    config
                        .residual_history_file
                        .as_ref()
                        .map(|p| format!(", output={}", p.display()))
                        .unwrap_or_default()
                );
            }
            println!("===============================================================");
            println!();
        }

        let base = Path::new(env!("CARGO_MANIFEST_DIR")).join("examples/mtx");
        let all_cases = [
            MatrixCase {
                file_name: "qc324.mtx",
                description: "qc324 (complex, general)",
                classification: "reference complex general case",
            },
            MatrixCase {
                file_name: "dwg961a.mtx",
                description: "dwg961a (complex, general)",
                classification: "stress/singular-or-bad-diagonal case",
            },
        ];
        let cases: Vec<MatrixCase> = match config.run_mode {
            RunMode::Correctness => all_cases.to_vec(),
            RunMode::Scalability => all_cases.to_vec(),
            RunMode::Robustness => all_cases.to_vec(),
        };

        if matches!(config.run_mode, RunMode::Correctness | RunMode::Scalability) {
            let (nx, ny) = config.poisson_grid;
            let generated_cases = [
                (
                    config.poisson_shift_real,
                    config.poisson_shift_imag,
                    "easy",
                    "generated complex shifted Poisson case",
                ),
                (
                    0.01,
                    0.25,
                    "hard",
                    "generated complex shifted Poisson hard-shift case",
                ),
            ];
            for (alpha, beta, difficulty, classification) in generated_cases {
                match build_shifted_poisson_case(
                    &comm,
                    nx,
                    ny,
                    alpha,
                    beta,
                    config.poisson_convection_x,
                    config.poisson_convection_y,
                ) {
                    Ok(problem) => execute_problem(
                        &format!(
                            "generated shifted Poisson {difficulty} (-Δ + ({alpha:.3} + {beta:.3}i)I, convection=({:.3},{:.3}))",
                            config.poisson_convection_x, config.poisson_convection_y
                        ),
                        classification,
                        problem,
                        &config,
                        rank,
                    ),
                    Err(err) if rank == 0 => {
                        println!(
                            "❌ Failed to build generated shifted Poisson {difficulty} case: {err}\n"
                        );
                    }
                    Err(_) => {}
                }
            }
        }

        for case in cases {
            let mat_name = case.file_name;
            let descr = case.description;
            let mat_path = base.join(mat_name);
            let available = mat_path.exists();
            if !available {
                if rank == 0 {
                    println!("⚠ Missing file {mat_name} for {descr}, skipping.\n");
                }
                continue;
            }

            let problem = match load_problem_complex(
                &mat_path,
                &comm,
                config.run_mode,
                config.mark_replicated_check,
            ) {
                Ok(p) => p,
                Err(err) => {
                    if rank == 0 {
                        println!("❌ Failed to load {descr}: {err}\n");
                    }
                    continue;
                }
            };

            execute_problem(descr, case.classification, problem, &config, rank);
        }

        if rank == 0 {
            println!("Example complete.");
            println!(
                "Final backend summary: {}.",
                if is_parallel {
                    "MPI run (see per-case backend labels for replicated/distributed mode)"
                } else {
                    "serial CSR"
                }
            );
        }

        Ok(())
    }

    fn execute_problem(
        descr: &str,
        classification: &str,
        problem: Problem,
        config: &BenchmarkConfig,
        rank: usize,
    ) {
        let rhs_norm2_local: f64 = problem.rhs.iter().map(|v| v.abs2()).sum();
        let rhs_norm = problem.comm.all_reduce_f64(rhs_norm2_local).sqrt();
        let runs = RunSpec::build_default_matrix(config, &problem);
        if rank == 0 {
            println!(
                "=== [{} mode] {descr} — {} ===",
                config.run_mode.label(),
                problem.backend_descr
            );
            println!("Classification: {classification}");
            println!(
                "Run backend: {} ({})",
                problem.backend.run_label(),
                problem.backend.details()
            );
            println!(
                "Benchmark/export label: {}",
                problem.backend.benchmark_export_label()
            );
            println!("Global DOFs: {}", problem.global_n);
            println!("Local DOFs (rank {rank}): {}", problem.local_n);
            println!(
                "Local row range (rank {rank}): [{}..{})",
                problem.global_row_start,
                problem.global_row_start + problem.local_n
            );
            if problem.comm.size() > 1 && problem.local_n == problem.global_n {
                println!("Note: replicated execution: MPI ranks are not sharing SpMV rows.");
            }
            if config.run_mode == RunMode::Correctness {
                println!(
                    "MPI scalability note: use distributed mode for scaling claims; correctness mode emphasizes solver validation."
                );
            }
            println!("‖b_unscaled‖₂ = {:.3e}", rhs_norm);
            println!(
                "{}",
                format_dist_policy_report(config, runs.first(), &problem)
            );
            println!("Residual semantics: rec/reported = solver recurrence/monitor residual.");
            if config.run_mode == RunMode::Correctness {
                println!(
                    "Verdict threshold: rel residual and generated x_err checks pass when <= {:.3e} ({}× rtol).",
                    correctness_verdict_tolerance(config),
                    CORRECTNESS_VERDICT_RTOL_MULTIPLIER
                );
                println!(
                    "                    True(dist) = ||b_unscaled - problem.op x_unscaled||₂ on owned rows, then reduced across ranks."
                );
                println!(
                    "                    True(global) = ||b_global - problem.global_csr x_global||₂ after gathering owned vectors in rank order."
                );
                println!(
                    "                    *(rel) columns divide by the corresponding unscaled RHS norm (stable under --row-scale)."
                );
            }
            println!(
                "FGMRES side policy: requested left/symmetric are normalized to effective right preconditioning."
            );
            if problem.rhs_source == RhsSource::GeneratedAOnes {
                println!(
                    "Generated A*ones reference: {}",
                    if problem.solution_reference.valid_for_x_error() {
                        "valid for x_err(rel) reporting"
                    } else {
                        "x_err(rel) suppressed; matrix has zero/missing rows or singularity indicators, so the solution may be nonunique"
                    }
                );
                if !problem.solution_reference.valid_for_x_error() {
                    println!("    → {}", problem.solution_reference.suppression_reason());
                }
            }
            if config.run_mode == RunMode::Correctness && config.mark_replicated_check {
                println!(
                    "Replicated check marker: ENABLED (metadata-only marker for cross-run comparison)."
                );
            }
            let include_dof_col = matches!(
                problem.backend,
                CsrBackend::Serial | CsrBackend::Distributed
            );
            if config.run_mode == RunMode::Scalability {
                println!(
                    "{:<7} {:<8} {:<6} {:<42} {:<36} {:<34} {:>9} {:>9} {:>7} {:>5} {:>5} {:>5} {:>17} {:>14} {:>12}",
                    "Op",
                    "Exec",
                    "PCdom",
                    "PC apply",
                    "Method",
                    "Effective policy",
                    "Med(s)",
                    "Min(s)",
                    "Iters",
                    "Rst",
                    "Inn",
                    "Pfb",
                    "Reds(obs/wait/model)",
                    "Rec/Reported",
                    "DOF/s"
                );
                println!("{}", "-".repeat(240));
            } else if include_dof_col {
                println!(
                    "{:<7} {:<8} {:<6} {:<42} {:<36} {:<34} {:<34} {:>9} {:>9} {:>9} {:>7} {:>5} {:>5} {:>5} {:>14} {:>14} {:>14} {:>14} {:>14} {:>14} {:>14} {:>8} {:>8} {:>8} {:>26} {:>12}",
                    "Op",
                    "Exec",
                    "PCdom",
                    "PC apply",
                    "Method",
                    "Requested policy",
                    "Effective policy",
                    "Setup(s)",
                    "Med(s)",
                    "Min(s)",
                    "Iters",
                    "Rst",
                    "Inn",
                    "Pfb",
                    "Rec/Reported",
                    "True(dist)",
                    "TrueDist(rel)",
                    "x_err(dist)",
                    "True(global)",
                    "TrueGlob(rel)",
                    "x_err(global)",
                    "DistOK",
                    "GlobalOK",
                    "XOK",
                    "Reason",
                    "DOF/s"
                );
                println!("{}", "-".repeat(408));
            } else {
                println!(
                    "{:<7} {:<8} {:<6} {:<42} {:<36} {:<34} {:<34} {:>9} {:>9} {:>9} {:>7} {:>5} {:>5} {:>5} {:>14} {:>14} {:>14} {:>14} {:>14} {:>14} {:>14} {:>8} {:>8} {:>8} {:>26}",
                    "Op",
                    "Exec",
                    "PCdom",
                    "PC apply",
                    "Method",
                    "Requested policy",
                    "Effective policy",
                    "Setup(s)",
                    "Med(s)",
                    "Min(s)",
                    "Iters",
                    "Rst",
                    "Inn",
                    "Pfb",
                    "Rec/Reported",
                    "True(dist)",
                    "TrueDist(rel)",
                    "x_err(dist)",
                    "True(global)",
                    "TrueGlob(rel)",
                    "x_err(global)",
                    "DistOK",
                    "GlobalOK",
                    "XOK",
                    "Reason"
                );
                println!("{}", "-".repeat(392));
            }
            println!(
                "Legend: Op=operator storage (csr-cx=complex CSR), Exec=execution backend (ser=serial, mpi-row=MPI row partition, mpi-repl=MPI replicated operator), PCdom=PC domain (full=global/full serial matrix ILU, own0=MPI owned block, overlap=0 local ILU/ASM, n/a=no ILU domain), PC apply=how the preconditioner is applied."
            );
        }

        for spec in runs {
            match run_once(&problem, &spec, config) {
                Ok(row) => {
                    if rank == 0 {
                        println!("{}", render_result_row(&row, config.run_mode, &problem));
                    }
                }
                Err(err) => {
                    if rank == 0 {
                        println!(
                            "{}",
                            render_failure_result_row(&spec, config.run_mode, &problem)
                        );
                        println!("    → {err}");
                    }
                }
            }
            problem.comm.barrier();
        }

        if rank == 0 {
            println!("{}", "=".repeat(96));
            println!();
        }
        problem.comm.barrier();
    }

    #[derive(Clone, Copy, Debug)]
    struct MatrixCase {
        file_name: &'static str,
        description: &'static str,
        classification: &'static str,
    }

    struct Problem {
        op: Arc<dyn KLinOp<Scalar = S>>,
        dist_plan_diagnostics: DistributedPlanDiagnostics,
        rhs: Vec<S>,
        rhs_source: RhsSource,
        solution_reference: SolutionReferenceDiagnostics,
        csr_for_pc: Arc<SparseCsrMatrix<S>>,
        global_csr: Arc<SparseCsrMatrix<S>>,
        local_rows_nnz: usize,
        zero_global_rows_local: usize,
        local_n: usize,
        global_n: usize,
        global_row_start: usize,
        comm: UniverseComm,
        backend: CsrBackend,
        backend_descr: String,
        generated_case: bool,
    }

    #[derive(Clone, Copy, Debug, Eq, PartialEq)]
    enum RhsSource {
        LoadedFromFile,
        GeneratedAOnes,
    }

    #[derive(Clone, Copy, Debug, Default, Eq, PartialEq)]
    struct SolutionReferenceDiagnostics {
        global_rows: usize,
        global_cols: usize,
        zero_row_count: usize,
        missing_diagonal_count: usize,
        tiny_diagonal_count: usize,
    }

    impl SolutionReferenceDiagnostics {
        fn valid_for_x_error(self) -> bool {
            self.global_rows == self.global_cols
                && self.zero_row_count == 0
                && self.missing_diagonal_count == 0
                && self.tiny_diagonal_count == 0
        }

        fn suppression_reason(self) -> String {
            let mut reasons = Vec::new();
            if self.global_rows != self.global_cols {
                reasons.push(format!(
                    "matrix is rectangular ({}x{})",
                    self.global_rows, self.global_cols
                ));
            }
            if self.zero_row_count > 0 {
                reasons.push(format!("{} zero/missing rows", self.zero_row_count));
            }
            if self.missing_diagonal_count > 0 {
                reasons.push(format!(
                    "{} structurally missing diagonals",
                    self.missing_diagonal_count
                ));
            }
            if self.tiny_diagonal_count > 0 {
                reasons.push(format!("{} tiny diagonals", self.tiny_diagonal_count));
            }
            if reasons.is_empty() {
                "generated A*ones reference is structurally valid".to_string()
            } else {
                format!(
                    "N/A (generated A*ones but matrix has {}; solution may be nonunique)",
                    reasons.join(", ")
                )
            }
        }
    }

    #[derive(Clone, Copy, Debug, Eq, PartialEq)]
    enum CsrBackend {
        Serial,
        Replicated,
        Distributed,
    }

    impl CsrBackend {
        fn run_label(self) -> &'static str {
            match self {
                Self::Serial => "serial CSR",
                Self::Replicated => "replicated CSR",
                Self::Distributed => "distributed CSR",
            }
        }

        fn details(self) -> &'static str {
            match self {
                Self::Serial => "single rank operator",
                Self::Replicated => "all ranks hold full matrix, identical solve",
                Self::Distributed => "row-partitioned operator",
            }
        }

        fn benchmark_export_label(self) -> &'static str {
            match self {
                Self::Serial => "serial",
                Self::Replicated => "replicated",
                Self::Distributed => "distributed",
            }
        }
    }

    const CORRECTNESS_VERDICT_RTOL_MULTIPLIER: f64 = 10.0;

    fn correctness_verdict_tolerance(bench_cfg: &BenchmarkConfig) -> f64 {
        bench_cfg.rtol * CORRECTNESS_VERDICT_RTOL_MULTIPLIER
    }

    struct ResultRow {
        operator_storage: &'static str,
        execution_backend: &'static str,
        pc_domain: &'static str,
        pc_apply: &'static str,
        method: String,
        requested_policy: String,
        effective_policy: String,
        setup_secs: f64,
        median_solve_secs: f64,
        min_solve_secs: f64,
        iterations: usize,
        reductions: usize,
        overlapped_reduction_waits: usize,
        model_predicted_reductions: Option<usize>,
        restart_count: Option<usize>,
        inner_iterations_last_cycle: Option<usize>,
        pipeline_fallbacks: Option<usize>,
        reported_residual: R,
        explicit_true_residual: Option<R>,
        explicit_true_residual_rel: Option<R>,
        x_error_rel: Option<R>,
        global_true_residual: Option<R>,
        global_true_residual_rel: Option<R>,
        global_x_error_rel: Option<R>,
        dist_ok: Option<bool>,
        global_ok: Option<bool>,
        x_ok: Option<bool>,
        reason: ConvergedReason,
        dof_per_sec: Option<f64>,
    }

    struct RunSpec {
        restart: usize,
        ksp: KspKind,
        variant: FgmresVariant,
        residual_check_policy: ResidualCheckPolicy,
        orthog: OrthogMethod,
        reorth: ReorthPolicy,
        pc_side: PcSide,
        pc: PcKind,
    }

    #[derive(Clone, Debug)]
    struct CsrForPcDiagnostics {
        nnz_local_block: usize,
        nnz_local_rows: usize,
        nnz_ratio: f64,
        zero_local_rows: usize,
        zero_global_rows: usize,
        structural_diag_missing: usize,
        diag_min_abs: f64,
        diag_max_abs: f64,
        diag_tiny_or_missing_count: usize,
    }

    #[derive(Clone, Debug)]
    struct RankSpread {
        global_sum: f64,
        rank_min: f64,
        rank_max: f64,
        rank0_local: f64,
    }

    #[derive(Clone, Debug)]
    struct CsrForPcDiagnosticsGlobal {
        nnz_local_block: RankSpread,
        nnz_local_rows: RankSpread,
        nnz_ratio: RankSpread,
        zero_local_rows_rank_count: f64,
        zero_global_rows: RankSpread,
        structural_diag_missing: RankSpread,
        diag_min_abs: RankSpread,
        diag_max_abs: RankSpread,
        diag_tiny_or_missing_count: RankSpread,
        pivot_perturbation_count: Option<RankSpread>,
    }

    #[derive(Clone, Copy, Debug, Eq, PartialEq)]
    enum PcKind {
        None,
        JacobiWeak,
        Ilu0Local,
        IlutLocal,
        MpiBlockJacobiIlu0Local,
        AsmIlu0Overlap { overlap: usize },
        RasIlu0Overlap { overlap: usize },
        RasIlukOverlap { k: usize, overlap: usize },
        LocalIluk { k: usize },
        ReplicatedFullIlu0,
        ReplicatedFullIluk { k: usize },
    }

    #[derive(Clone, Copy, Debug, Eq, PartialEq)]
    enum PcDispatchBranch {
        None,
        JacobiWeak,
        Ilu0Local,
        ReplicatedFullIlu,
        IlutLocal,
        OverlapIlu,
    }

    #[derive(Clone, Debug)]
    struct BenchmarkConfig {
        run_mode: RunMode,
        warmup_runs: usize,
        measured_runs: usize,
        rtol: f64,
        atol: f64,
        maxits: usize,
        restarts: Vec<usize>,
        ksp: KspKind,
        pcs: Vec<PcKind>,
        include_restart_200: bool,
        include_ilut_real_projection_fallback: bool,
        include_overlap_pcs: bool,
        variants: Vec<FgmresVariant>,
        orthogs: Vec<OrthogMethod>,
        reorths: Vec<ReorthPolicy>,
        dist_policy: DistPolicyMode,
        allow_stagnation_fallback: bool,
        min_inner_before_fallback: usize,
        mark_replicated_check: bool,
        residual_history: bool,
        residual_history_file: Option<PathBuf>,
        residual_history_force: bool,
        fgmres_haptol: f64,
        row_scale: bool,
        row_scale_tiny: f64,
        ilu_reordering: ReorderingOptions,
        poisson_grid: (usize, usize),
        poisson_shift_real: f64,
        poisson_shift_imag: f64,
        poisson_convection_x: f64,
        poisson_convection_y: f64,
    }

    #[derive(Clone, Copy, Debug, Eq, PartialEq)]
    enum RunMode {
        Correctness,
        Scalability,
        Robustness,
    }

    impl RunMode {
        fn default_dist_policy(self) -> DistPolicyMode {
            match self {
                Self::Correctness => DistPolicyMode::Off,
                Self::Scalability => DistPolicyMode::Auto,
                // Robustness keeps the requested restart/variant unless explicitly opted in.
                Self::Robustness => DistPolicyMode::Off,
            }
        }

        fn label(self) -> &'static str {
            match self {
                Self::Correctness => "correctness",
                Self::Scalability => "scalability",
                Self::Robustness => "robustness",
            }
        }

        fn parse(token: &str) -> Result<Self, KError> {
            match token.trim().to_ascii_lowercase().as_str() {
                "correctness" => Ok(Self::Correctness),
                "scalability" => Ok(Self::Scalability),
                "robustness" => Ok(Self::Robustness),
                other => Err(KError::InvalidInput(format!(
                    "invalid run mode '{other}', expected correctness|scalability|robustness"
                ))),
            }
        }
    }

    #[derive(Clone, Copy, Debug, Eq, PartialEq)]
    enum DistPolicyMode {
        Off,
        Auto,
    }

    impl DistPolicyMode {
        fn label(self) -> &'static str {
            match self {
                Self::Off => "off",
                Self::Auto => "auto",
            }
        }

        fn parse(token: &str) -> Result<Self, KError> {
            match token.trim().to_ascii_lowercase().as_str() {
                "off" => Ok(Self::Off),
                "auto" => Ok(Self::Auto),
                other => Err(KError::InvalidInput(format!(
                    "invalid DistCSR policy '{other}', expected off|auto"
                ))),
            }
        }
    }

    #[derive(Clone, Copy, Debug, Eq, PartialEq)]
    enum KspKind {
        Gmres,
        Fgmres,
    }

    impl KspKind {
        fn label(self) -> &'static str {
            match self {
                Self::Gmres => "GMRES",
                Self::Fgmres => "FGMRES",
            }
        }

        fn parse(token: &str) -> Result<Self, KError> {
            match token.trim().to_ascii_lowercase().as_str() {
                "gmres" => Ok(Self::Gmres),
                "fgmres" => Ok(Self::Fgmres),
                other => Err(KError::InvalidInput(format!(
                    "invalid KSP '{other}', expected gmres|fgmres"
                ))),
            }
        }

        fn for_pc(self, pc: PcKind) -> Self {
            if pc.requires_flexible_right_preconditioning() {
                Self::Fgmres
            } else {
                self
            }
        }
    }

    impl Default for BenchmarkConfig {
        fn default() -> Self {
            Self {
                run_mode: RunMode::Correctness,
                warmup_runs: 1,
                measured_runs: 5,
                rtol: 1e-8,
                atol: 1e-12,
                maxits: 500,
                restarts: vec![50, 100, 150],
                ksp: KspKind::Gmres,
                pcs: Vec::new(),
                include_restart_200: false,
                include_ilut_real_projection_fallback: false,
                include_overlap_pcs: false,
                variants: vec![FgmresVariant::Classical],
                orthogs: vec![OrthogMethod::ClassicalGS],
                reorths: vec![ReorthPolicy::IfNeeded],
                dist_policy: RunMode::Correctness.default_dist_policy(),
                allow_stagnation_fallback: false,
                min_inner_before_fallback: 8,
                mark_replicated_check: false,
                residual_history: false,
                residual_history_file: None,
                residual_history_force: false,
                fgmres_haptol: 1e-30,
                row_scale: false,
                row_scale_tiny: 1e-15,
                ilu_reordering: ReorderingOptions::default(),
                poisson_grid: (16, 16),
                poisson_shift_real: 1.0,
                poisson_shift_imag: 0.25,
                poisson_convection_x: 0.0,
                poisson_convection_y: 0.0,
            }
        }
    }

    impl BenchmarkConfig {
        fn from_env_args() -> Result<Self, KError> {
            Self::from_args(env::args().skip(1))
        }

        fn from_args<I>(args: I) -> Result<Self, KError>
        where
            I: IntoIterator<Item = String>,
        {
            let mut cfg = Self::default();
            let mut dist_policy_explicit = false;
            let mut args = args.into_iter().peekable();
            while let Some(arg) = args.next() {
                match arg.as_str() {
                    "--warmup-runs" => {
                        let Some(v) = args.next() else {
                            return Err(KError::InvalidInput(
                                "missing value for --warmup-runs".into(),
                            ));
                        };
                        cfg.warmup_runs = parse_positive_usize("--warmup-runs", &v)?;
                    }
                    "--measured-runs" => {
                        let Some(v) = args.next() else {
                            return Err(KError::InvalidInput(
                                "missing value for --measured-runs".into(),
                            ));
                        };
                        cfg.measured_runs = parse_positive_usize("--measured-runs", &v)?;
                    }
                    "--help" | "-h" => {
                        if cfg!(feature = "mpi") {
                            println!(
                                "Usage: cargo mpirun -n <ranks> --example complex_matrix_market_demo --features complex,mpi,mpi_examples -- [--mode correctness|scalability|robustness] [--ksp gmres|fgmres] [--dist-policy off|auto] [--mark-replicated-check] [--warmup-runs N] [--measured-runs N] [--rtol F] [--atol F] [--maxits N] [--restarts csv] [--pcs csv] [--include-restart-200] [--include-ilut-real-projection-fallback] [--include-overlap-pcs] [--allow-stagnation-fallback] [--min-inner-before-fallback N] [--fgmres-variant csv] [--fgmres-orthog csv] [--fgmres-reorth csv] [--fgmres-haptol F] [--residual-history] [--residual-history-file <path>] [--residual-history-force] [--row-scale [tiny]] [--ilu-reordering none|rcm|amd[:nonsym]] [--poisson-grid NXxNY] [--poisson-shift-real F] [--poisson-shift-imag F] [--poisson-convection CX,CY]
Defaults: --ksp gmres for fixed linear PCs; flexible/nonlinear PC modes select FGMRES. --dist-policy is off for correctness and robustness, auto for scalability.
--include-overlap-pcs is active only for MPI correctness sweeps; it keeps block-jacobi-ilu0-overlap0 and appends asm-ilu0-overlap1, ras-ilu0-overlap1, ras-iluk1-overlap1 so GlobalOK/XOK verdicts are visible before performance comparisons."
                            );
                        } else {
                            println!(
                                "Usage: cargo run --example complex_matrix_market_demo --features complex -- [--mode correctness|scalability|robustness] [--ksp gmres|fgmres] [--dist-policy off|auto] [--mark-replicated-check] [--warmup-runs N] [--measured-runs N] [--rtol F] [--atol F] [--maxits N] [--restarts csv] [--pcs csv] [--include-restart-200] [--include-ilut-real-projection-fallback] [--include-overlap-pcs] [--allow-stagnation-fallback] [--min-inner-before-fallback N] [--fgmres-variant csv] [--fgmres-orthog csv] [--fgmres-reorth csv] [--fgmres-haptol F] [--residual-history] [--residual-history-file <path>] [--residual-history-force] [--row-scale [tiny]] [--ilu-reordering none|rcm|amd[:nonsym]] [--poisson-grid NXxNY] [--poisson-shift-real F] [--poisson-shift-imag F] [--poisson-convection CX,CY]
Defaults: --ksp gmres for fixed linear PCs; flexible/nonlinear PC modes select FGMRES. --dist-policy is off for correctness and robustness, auto for scalability.
--include-overlap-pcs is active only for MPI correctness sweeps; it keeps block-jacobi-ilu0-overlap0 and appends asm-ilu0-overlap1, ras-ilu0-overlap1, ras-iluk1-overlap1 so GlobalOK/XOK verdicts are visible before performance comparisons."
                            );
                        }
                        std::process::exit(0);
                    }
                    "--mode" => {
                        let Some(v) = args.next() else {
                            return Err(KError::InvalidInput("missing value for --mode".into()));
                        };
                        cfg.run_mode = RunMode::parse(&v)?;
                    }
                    "--ksp" => {
                        let Some(v) = args.next() else {
                            return Err(KError::InvalidInput("missing value for --ksp".into()));
                        };
                        cfg.ksp = KspKind::parse(&v)?;
                    }
                    "--dist-policy" => {
                        let Some(v) = args.next() else {
                            return Err(KError::InvalidInput(
                                "missing value for --dist-policy".into(),
                            ));
                        };
                        cfg.dist_policy = DistPolicyMode::parse(&v)?;
                        dist_policy_explicit = true;
                    }
                    "--poisson-grid" => {
                        let Some(v) = args.next() else {
                            return Err(KError::InvalidInput(
                                "missing value for --poisson-grid".into(),
                            ));
                        };
                        cfg.poisson_grid = parse_grid_dims("--poisson-grid", &v)?;
                    }
                    "--poisson-shift-real" => {
                        let Some(v) = args.next() else {
                            return Err(KError::InvalidInput(
                                "missing value for --poisson-shift-real".into(),
                            ));
                        };
                        cfg.poisson_shift_real = parse_finite_f64("--poisson-shift-real", &v)?;
                    }
                    "--poisson-shift-imag" => {
                        let Some(v) = args.next() else {
                            return Err(KError::InvalidInput(
                                "missing value for --poisson-shift-imag".into(),
                            ));
                        };
                        cfg.poisson_shift_imag = parse_finite_f64("--poisson-shift-imag", &v)?;
                    }
                    "--poisson-convection" => {
                        let Some(v) = args.next() else {
                            return Err(KError::InvalidInput(
                                "missing value for --poisson-convection".into(),
                            ));
                        };
                        (cfg.poisson_convection_x, cfg.poisson_convection_y) =
                            parse_f64_pair("--poisson-convection", &v)?;
                    }
                    "--mark-replicated-check" => {
                        cfg.mark_replicated_check = true;
                    }
                    "--correctness-replicated-check" => {
                        return Err(KError::InvalidInput(
                            "--correctness-replicated-check is metadata-only and has been renamed to --mark-replicated-check".into(),
                        ));
                    }
                    "--restarts" => {
                        let Some(v) = args.next() else {
                            return Err(KError::InvalidInput(
                                "missing value for --restarts".into(),
                            ));
                        };
                        cfg.restarts = parse_usize_csv("--restarts", &v)?;
                    }
                    "--rtol" => {
                        let Some(v) = args.next() else {
                            return Err(KError::InvalidInput("missing value for --rtol".into()));
                        };
                        cfg.rtol = parse_positive_f64("--rtol", &v)?;
                    }
                    "--atol" => {
                        let Some(v) = args.next() else {
                            return Err(KError::InvalidInput("missing value for --atol".into()));
                        };
                        cfg.atol = parse_positive_f64("--atol", &v)?;
                    }
                    "--maxits" => {
                        let Some(v) = args.next() else {
                            return Err(KError::InvalidInput("missing value for --maxits".into()));
                        };
                        cfg.maxits = parse_positive_usize("--maxits", &v)?;
                    }
                    "--pcs" => {
                        let Some(v) = args.next() else {
                            return Err(KError::InvalidInput("missing value for --pcs".into()));
                        };
                        cfg.pcs = parse_pc_csv("--pcs", &v)?;
                    }
                    "--include-restart-200" => {
                        cfg.include_restart_200 = true;
                    }
                    "--include-ilut-real-projection-fallback" => {
                        cfg.include_ilut_real_projection_fallback = true;
                    }
                    "--include-overlap-pcs" => {
                        cfg.include_overlap_pcs = true;
                    }
                    "--ilu-reordering" | "--pc-ilu-reordering-type" => {
                        let Some(v) = args.next() else {
                            return Err(KError::InvalidInput(
                                "missing value for --ilu-reordering".into(),
                            ));
                        };
                        cfg.ilu_reordering = parse_ilu_reordering(&v)?;
                    }
                    "--allow-stagnation-fallback" => {
                        cfg.allow_stagnation_fallback = true;
                    }
                    "--min-inner-before-fallback" => {
                        let Some(v) = args.next() else {
                            return Err(KError::InvalidInput(
                                "missing value for --min-inner-before-fallback".into(),
                            ));
                        };
                        cfg.min_inner_before_fallback =
                            parse_positive_usize("--min-inner-before-fallback", &v)?;
                    }
                    "--fgmres-variant" => {
                        let Some(v) = args.next() else {
                            return Err(KError::InvalidInput(
                                "missing value for --fgmres-variant".into(),
                            ));
                        };
                        cfg.variants = parse_variant_csv("--fgmres-variant", &v)?;
                    }
                    "--fgmres-orthog" => {
                        let Some(v) = args.next() else {
                            return Err(KError::InvalidInput(
                                "missing value for --fgmres-orthog".into(),
                            ));
                        };
                        cfg.orthogs = parse_orthog_csv("--fgmres-orthog", &v)?;
                    }
                    "--residual-history" => {
                        cfg.residual_history = true;
                    }
                    "--residual-history-file" => {
                        let Some(v) = args.next() else {
                            return Err(KError::InvalidInput(
                                "missing value for --residual-history-file".into(),
                            ));
                        };
                        cfg.residual_history_file = Some(PathBuf::from(v));
                    }
                    "--residual-history-force" => {
                        cfg.residual_history_force = true;
                    }
                    "--row-scale" => {
                        cfg.row_scale = true;
                        if let Some(next) = args.peek() {
                            if !next.starts_with("--") {
                                let tiny = args.next().ok_or_else(|| {
                                    KError::InvalidInput(
                                        "failed to read --row-scale optional tiny".into(),
                                    )
                                })?;
                                cfg.row_scale_tiny = parse_positive_f64("--row-scale", &tiny)?;
                            }
                        }
                    }
                    "--fgmres-reorth" => {
                        let Some(v) = args.next() else {
                            return Err(KError::InvalidInput(
                                "missing value for --fgmres-reorth".into(),
                            ));
                        };
                        cfg.reorths = parse_reorth_csv("--fgmres-reorth", &v)?;
                    }
                    "--fgmres-haptol" => {
                        let Some(v) = args.next() else {
                            return Err(KError::InvalidInput(
                                "missing value for --fgmres-haptol".into(),
                            ));
                        };
                        cfg.fgmres_haptol = parse_positive_finite_f64("--fgmres-haptol", &v)?;
                    }
                    _ => {
                        return Err(KError::InvalidInput(format!("unknown argument: {arg}")));
                    }
                }
            }
            if cfg.measured_runs == 0 {
                return Err(KError::InvalidInput(
                    "--measured-runs must be at least 1".into(),
                ));
            }
            if !dist_policy_explicit {
                cfg.dist_policy = cfg.run_mode.default_dist_policy();
            }
            if cfg.include_restart_200 && !cfg.restarts.contains(&200) {
                cfg.restarts.push(200);
            }
            if cfg.run_mode == RunMode::Scalability {
                cfg.mark_replicated_check = false;
            }
            if cfg.pcs.iter().any(|pc| pc.is_ilut()) {
                eprintln!(
                    "⚠ --pcs includes ILUT for a complex run: current ILUT path is a degraded real projection and is not trusted for complex robustness benchmarking."
                );
            }
            if cfg.include_ilut_real_projection_fallback && !cfg.pcs.iter().any(|pc| pc.is_ilut()) {
                eprintln!(
                    "⚠ Including {} in the default complex benchmark set: current ILUT path is a degraded real projection and is not trusted for complex robustness benchmarking.",
                    ILUT_REAL_PROJECTION_FALLBACK_LABEL
                );
            }
            if cfg.residual_history
                && cfg.run_mode != RunMode::Correctness
                && !cfg.residual_history_force
            {
                cfg.residual_history = false;
            }
            Ok(cfg)
        }
    }

    fn parse_grid_dims(flag: &str, value: &str) -> Result<(usize, usize), KError> {
        let Some((nx, ny)) = value.split_once(['x', 'X']) else {
            return Err(KError::InvalidInput(format!(
                "invalid value '{value}' for {flag}, expected NXxNY"
            )));
        };
        let nx = parse_positive_usize(flag, nx)?;
        let ny = parse_positive_usize(flag, ny)?;
        if nx == 0 || ny == 0 {
            return Err(KError::InvalidInput(format!(
                "invalid value '{value}' for {flag}, dimensions must be positive"
            )));
        }
        Ok((nx, ny))
    }

    fn parse_positive_usize(flag: &str, value: &str) -> Result<usize, KError> {
        value.parse::<usize>().map_err(|_| {
            KError::InvalidInput(format!(
                "invalid value '{value}' for {flag}, expected non-negative integer"
            ))
        })
    }

    fn parse_positive_f64(flag: &str, value: &str) -> Result<f64, KError> {
        let val = value.parse::<f64>().map_err(|_| {
            KError::InvalidInput(format!(
                "invalid value '{value}' for {flag}, expected non-negative float"
            ))
        })?;
        if val < 0.0 {
            return Err(KError::InvalidInput(format!(
                "invalid value '{value}' for {flag}, expected non-negative float"
            )));
        }
        Ok(val)
    }

    fn parse_finite_f64(flag: &str, value: &str) -> Result<f64, KError> {
        let val = value.parse::<f64>().map_err(|_| {
            KError::InvalidInput(format!(
                "invalid value '{value}' for {flag}, expected finite float"
            ))
        })?;
        if !val.is_finite() {
            return Err(KError::InvalidInput(format!(
                "invalid value '{value}' for {flag}, expected finite float"
            )));
        }
        Ok(val)
    }

    fn parse_f64_pair(flag: &str, value: &str) -> Result<(f64, f64), KError> {
        let Some((x, y)) = value.split_once(',') else {
            return Err(KError::InvalidInput(format!(
                "invalid value '{value}' for {flag}, expected X,Y"
            )));
        };
        let x = parse_finite_f64(flag, x.trim())?;
        let y = parse_finite_f64(flag, y.trim())?;
        Ok((x, y))
    }

    fn parse_positive_finite_f64(flag: &str, value: &str) -> Result<f64, KError> {
        let val = value.parse::<f64>().map_err(|_| {
            KError::InvalidInput(format!(
                "invalid value '{value}' for {flag}, expected positive finite float"
            ))
        })?;
        if !val.is_finite() || val <= 0.0 {
            return Err(KError::InvalidInput(format!(
                "invalid value '{value}' for {flag}, expected positive finite float"
            )));
        }
        Ok(val)
    }

    fn parse_usize_csv(flag: &str, value: &str) -> Result<Vec<usize>, KError> {
        let vals = value
            .split(',')
            .map(str::trim)
            .filter(|s| !s.is_empty())
            .map(|v| parse_positive_usize(flag, v))
            .collect::<Result<Vec<_>, _>>()?;
        if vals.is_empty() {
            return Err(KError::InvalidInput(format!(
                "{flag} expects at least one integer value"
            )));
        }
        Ok(vals)
    }

    fn parse_variant(token: &str) -> Result<FgmresVariant, KError> {
        match token.trim().to_ascii_lowercase().as_str() {
            "classical" => Ok(FgmresVariant::Classical),
            "pipelined" => Ok(FgmresVariant::Pipelined),
            other => Err(KError::InvalidInput(format!(
                "invalid fgmres variant '{other}', expected classical|pipelined"
            ))),
        }
    }

    fn parse_orthog(token: &str) -> Result<OrthogMethod, KError> {
        match token.trim().to_ascii_lowercase().as_str() {
            "classical" | "cgs" | "cgs_refined" | "cgs-refined" | "refined" => {
                Ok(OrthogMethod::ClassicalGS)
            }
            "mgs" | "modified" => Ok(OrthogMethod::ModifiedGS),
            other => Err(KError::InvalidInput(format!(
                "invalid orthog '{other}', expected classical|cgs|cgs_refined|mgs|modified"
            ))),
        }
    }

    fn parse_reorth(token: &str) -> Result<ReorthPolicy, KError> {
        match token.trim().to_ascii_lowercase().as_str() {
            "never" => Ok(ReorthPolicy::Never),
            "ifneeded" | "if-needed" => Ok(ReorthPolicy::IfNeeded),
            "always" => Ok(ReorthPolicy::Always),
            other => Err(KError::InvalidInput(format!(
                "invalid reorth '{other}', expected never|ifneeded|always"
            ))),
        }
    }

    fn parse_variant_csv(flag: &str, value: &str) -> Result<Vec<FgmresVariant>, KError> {
        parse_csv(flag, value, parse_variant)
    }

    fn parse_orthog_csv(flag: &str, value: &str) -> Result<Vec<OrthogMethod>, KError> {
        parse_csv(flag, value, parse_orthog)
    }

    fn parse_reorth_csv(flag: &str, value: &str) -> Result<Vec<ReorthPolicy>, KError> {
        parse_csv(flag, value, parse_reorth)
    }

    fn parse_ilu_reordering(value: &str) -> Result<ReorderingOptions, KError> {
        let normalized = value.trim().to_ascii_lowercase();
        let (base, nonsym) = if let Some(base) = normalized.strip_suffix(":nonsym") {
            (base, true)
        } else if let Some(base) = normalized.strip_suffix("_nonsym") {
            (base, true)
        } else {
            (normalized.as_str(), false)
        };
        let kind = match base {
            "none" | "natural" => ReorderingKind::None,
            "rcm" => ReorderingKind::Rcm,
            "amd" => ReorderingKind::Amd,
            other => {
                return Err(KError::InvalidInput(format!(
                    "invalid ILU reordering '{other}', expected none|rcm|amd with optional :nonsym suffix"
                )));
            }
        };
        Ok(ReorderingOptions {
            kind,
            symmetric: !nonsym,
            deterministic: true,
        })
    }

    fn parse_pc(token: &str) -> Result<PcKind, KError> {
        let token_norm = token.trim().to_ascii_lowercase();
        if let Some(k_str) = token_norm.strip_prefix("local-iluk:") {
            return Ok(PcKind::LocalIluk {
                k: parse_positive_usize("local-iluk:k", k_str)?,
            });
        }
        if let Some(k_str) = token_norm.strip_prefix("block-jacobi-iluk:") {
            return Ok(PcKind::LocalIluk {
                k: parse_positive_usize("block-jacobi-iluk:k", k_str)?,
            });
        }
        if let Some(k_str) = token_norm.strip_prefix("block-jacobi-iluk") {
            let k_str = k_str.strip_suffix("-overlap0").unwrap_or(k_str);
            return Ok(PcKind::LocalIluk {
                k: parse_positive_usize("block-jacobi-iluk:k", k_str)?,
            });
        }
        if let Some(k_str) = token_norm.strip_prefix("replicated-iluk:") {
            return Ok(PcKind::ReplicatedFullIluk {
                k: parse_positive_usize("replicated-iluk:k", k_str)?,
            });
        }
        match token_norm.as_str() {
            "asm-ilu0-overlap1" => Ok(PcKind::AsmIlu0Overlap { overlap: 1 }),
            "asm-ilu0-overlap2" => Ok(PcKind::AsmIlu0Overlap { overlap: 2 }),
            "ras-ilu0-overlap1" => Ok(PcKind::RasIlu0Overlap { overlap: 1 }),
            "ras-iluk1-overlap1" => Ok(PcKind::RasIlukOverlap { k: 1, overlap: 1 }),
            "none" | "off" => Ok(PcKind::None),
            "jacobi" | "jacobi-weak" | "weak-jacobi" => Ok(PcKind::JacobiWeak),
            "ilu0" | "ilu0-local" | "local-ilu0" => Ok(PcKind::Ilu0Local),
            "ilut" | "ilut-local" | "local-ilut" => Ok(PcKind::IlutLocal),
            "replicated-ilu0" | "replicated-full-ilu0" => Ok(PcKind::ReplicatedFullIlu0),
            "mpi-block-jacobi-ilu0"
            | "block-jacobi-ilu0"
            | "block-jacobi-ilu0-overlap0"
            | "mpi-block-ilu0" => Ok(PcKind::MpiBlockJacobiIlu0Local),
            other => Err(KError::InvalidInput(format!(
                "invalid pc '{other}', expected none|jacobi|local-ilu0|local-ilut|local-iluk:<k>|block-jacobi-ilu0-overlap0|block-jacobi-iluk1|block-jacobi-iluk1-overlap0|block-jacobi-iluk:<k>|block-jacobi-iluk<k>-overlap0|replicated-ilu0|replicated-full-ilu0|replicated-iluk:<k>|mpi-block-jacobi-ilu0|asm-ilu0-overlap1|asm-ilu0-overlap2|ras-ilu0-overlap1|ras-iluk1-overlap1"
            ))),
        }
    }

    fn parse_pc_csv(flag: &str, value: &str) -> Result<Vec<PcKind>, KError> {
        parse_csv(flag, value, parse_pc)
    }

    const MPI_CORRECTNESS_OVERLAP_PCS: [PcKind; 3] = [
        PcKind::AsmIlu0Overlap { overlap: 1 },
        PcKind::RasIlu0Overlap { overlap: 1 },
        PcKind::RasIlukOverlap { k: 1, overlap: 1 },
    ];

    fn should_include_overlap_pcs(cfg: &BenchmarkConfig, comm_size: usize) -> bool {
        cfg.include_overlap_pcs && cfg.run_mode == RunMode::Correctness && comm_size > 1
    }

    fn push_pc_unique(pcs: &mut Vec<PcKind>, pc: PcKind) {
        if !pcs.contains(&pc) {
            pcs.push(pc);
        }
    }

    fn append_mpi_correctness_overlap_pcs(pcs: &mut Vec<PcKind>) {
        push_pc_unique(pcs, PcKind::MpiBlockJacobiIlu0Local);
        for pc in MPI_CORRECTNESS_OVERLAP_PCS {
            push_pc_unique(pcs, pc);
        }
    }

    fn parse_csv<T, F>(flag: &str, value: &str, mut parser: F) -> Result<Vec<T>, KError>
    where
        F: FnMut(&str) -> Result<T, KError>,
    {
        let vals = value
            .split(',')
            .map(str::trim)
            .filter(|s| !s.is_empty())
            .map(&mut parser)
            .collect::<Result<Vec<_>, _>>()?;
        if vals.is_empty() {
            return Err(KError::InvalidInput(format!(
                "{flag} expects at least one value"
            )));
        }
        Ok(vals)
    }

    impl RunSpec {
        fn build_default_matrix(cfg: &BenchmarkConfig, problem: &Problem) -> Vec<Self> {
            let mut pcs = if cfg.pcs.is_empty() && problem.generated_case {
                vec![
                    PcKind::None,
                    PcKind::JacobiWeak,
                    if problem.comm.size() > 1 {
                        PcKind::MpiBlockJacobiIlu0Local
                    } else {
                        PcKind::Ilu0Local
                    },
                    PcKind::LocalIluk { k: 1 },
                ]
            } else if cfg.pcs.is_empty() {
                match cfg.run_mode {
                    RunMode::Correctness => {
                        if problem.comm.size() > 1 {
                            vec![
                                PcKind::ReplicatedFullIlu0,
                                PcKind::None,
                                PcKind::JacobiWeak,
                                PcKind::MpiBlockJacobiIlu0Local,
                            ]
                        } else {
                            vec![PcKind::Ilu0Local, PcKind::JacobiWeak, PcKind::None]
                        }
                    }
                    RunMode::Scalability => {
                        if problem.comm.size() > 1 {
                            vec![PcKind::MpiBlockJacobiIlu0Local]
                        } else {
                            vec![PcKind::Ilu0Local]
                        }
                    }
                    RunMode::Robustness => {
                        if problem.comm.size() > 1 {
                            vec![PcKind::MpiBlockJacobiIlu0Local, PcKind::JacobiWeak]
                        } else {
                            vec![PcKind::Ilu0Local, PcKind::JacobiWeak]
                        }
                    }
                }
            } else {
                cfg.pcs.clone()
            };
            if cfg.pcs.is_empty() && cfg.include_ilut_real_projection_fallback {
                pcs.push(PcKind::IlutLocal);
            }
            if should_include_overlap_pcs(cfg, problem.comm.size()) {
                append_mpi_correctness_overlap_pcs(&mut pcs);
            }
            let mut runs = Vec::new();
            for &restart in &cfg.restarts {
                for &variant in &cfg.variants {
                    for &orthog in &cfg.orthogs {
                        for &reorth in &cfg.reorths {
                            for &pc in &pcs {
                                if pc == PcKind::MpiBlockJacobiIlu0Local && problem.comm.size() <= 1
                                {
                                    continue;
                                }
                                runs.push(Self {
                                    restart,
                                    ksp: cfg.ksp.for_pc(pc),
                                    variant,
                                    residual_check_policy: match cfg.run_mode {
                                        RunMode::Correctness => ResidualCheckPolicy::OnConvergence,
                                        RunMode::Scalability => ResidualCheckPolicy::RestartOnly,
                                        RunMode::Robustness => ResidualCheckPolicy::RestartOnly,
                                    },
                                    orthog,
                                    reorth,
                                    pc_side: PcSide::Right,
                                    pc,
                                });
                            }
                        }
                    }
                }
            }
            runs
        }

        fn method_label(&self) -> String {
            let effective_side = normalized_fgmres_side(self.pc_side);
            let side_desc = if effective_side == self.pc_side {
                format!("{}", pc_side_label(self.pc_side))
            } else {
                format!(
                    "{}→{} (normalized)",
                    pc_side_label(self.pc_side),
                    pc_side_label(effective_side)
                )
            };
            format!(
                "{}+{} [m={}, v={}, reschk={}, orth={}, reorth={}, pc=requested {}, effective {}]",
                self.ksp.label(),
                self.pc.label(),
                self.restart,
                variant_label(self.variant),
                residual_check_policy_label(self.residual_check_policy),
                orthog_label(self.orthog),
                reorth_label(self.reorth),
                pc_side_label(self.pc_side),
                side_desc,
            )
        }

        fn requested_policy_label(&self) -> String {
            format!(
                "ksp={}, variant={}, restart={}, residual-check={}",
                self.ksp.label(),
                variant_label(self.variant),
                self.restart,
                residual_check_policy_label(self.residual_check_policy)
            )
        }
    }

    fn operator_storage_label(_problem: &Problem) -> &'static str {
        "csr-cx"
    }

    fn execution_backend_label(problem: &Problem) -> &'static str {
        match problem.backend {
            CsrBackend::Serial => "ser",
            CsrBackend::Replicated => "mpi-repl",
            CsrBackend::Distributed => "mpi-row",
        }
    }

    fn pc_domain_label(pc: PcKind, problem: &Problem) -> &'static str {
        pc_domain_label_for_size(pc, problem.comm.size())
    }

    fn pc_domain_label_for_size(pc: PcKind, comm_size: usize) -> &'static str {
        match pc {
            PcKind::ReplicatedFullIlu0 | PcKind::ReplicatedFullIluk { .. } => "full",
            PcKind::Ilu0Local | PcKind::IlutLocal | PcKind::LocalIluk { .. } => {
                if comm_size <= 1 {
                    "full"
                } else {
                    "own0"
                }
            }
            PcKind::MpiBlockJacobiIlu0Local => "own0",
            PcKind::AsmIlu0Overlap { overlap }
            | PcKind::RasIlu0Overlap { overlap }
            | PcKind::RasIlukOverlap { overlap, .. } => match overlap {
                1 => "own+ghost1",
                2 => "own+ghost2",
                _ => "own+ghostN",
            },
            PcKind::None | PcKind::JacobiWeak => "n/a",
        }
    }

    fn pc_apply_label(pc: PcKind, problem: &Problem) -> &'static str {
        pc_apply_label_for_size(pc, problem.comm.size())
    }

    fn pc_apply_label_for_size(pc: PcKind, comm_size: usize) -> &'static str {
        match pc {
            PcKind::ReplicatedFullIlu0 | PcKind::ReplicatedFullIluk { .. } => {
                if comm_size > 1 {
                    "replicated/full/allgather, not scalable"
                } else {
                    "replicated/full"
                }
            }
            PcKind::Ilu0Local
            | PcKind::IlutLocal
            | PcKind::MpiBlockJacobiIlu0Local
            | PcKind::LocalIluk { .. } => "owned-block/overlap0",
            PcKind::AsmIlu0Overlap { .. } => "asm/overlap-gather/owned-output",
            PcKind::RasIlu0Overlap { .. } | PcKind::RasIlukOverlap { .. } => {
                "ras/overlap-gather/restricted-owned"
            }
            PcKind::JacobiWeak => "jacobi/local",
            PcKind::None => "none",
        }
    }

    const ILUT_REAL_PROJECTION_FALLBACK_LABEL: &str = "ILUT(real-projection fallback)";

    impl PcKind {
        fn label(&self) -> String {
            match self {
                Self::None => "none (unpreconditioned reference)".to_string(),
                Self::JacobiWeak => "jacobi (weak baseline)".to_string(),
                Self::Ilu0Local => "block-jacobi-ilu0-overlap0".to_string(),
                Self::IlutLocal => {
                    "local ILUT(real-projection fallback) [degraded/provisional complex path: real-projection fallback; not trusted for complex robustness benchmarking]".to_string()
                }
                Self::MpiBlockJacobiIlu0Local => {
                    "block-jacobi-ilu0-overlap0 [mpi spelling alias]".to_string()
                }
                Self::AsmIlu0Overlap { overlap } => format!("asm-ilu0-overlap{overlap}"),
                Self::RasIlu0Overlap { overlap } => format!("ras-ilu0-overlap{overlap}"),
                Self::RasIlukOverlap { k, overlap } => format!("ras-iluk{k}-overlap{overlap}"),
                Self::LocalIluk { k } => format!("block-jacobi-iluk{k}-overlap0"),
                Self::ReplicatedFullIlu0 => {
                    "replicated full ILU(0) [correctness only, not scalable]".to_string()
                }
                Self::ReplicatedFullIluk { .. } => {
                    "replicated full ILU(k) [correctness only, not scalable]".to_string()
                }
            }
        }

        fn dispatch_branch(self) -> PcDispatchBranch {
            match self {
                Self::None => PcDispatchBranch::None,
                Self::JacobiWeak => PcDispatchBranch::JacobiWeak,
                Self::Ilu0Local | Self::MpiBlockJacobiIlu0Local | Self::LocalIluk { .. } => {
                    PcDispatchBranch::Ilu0Local
                }
                Self::ReplicatedFullIlu0 | Self::ReplicatedFullIluk { .. } => {
                    PcDispatchBranch::ReplicatedFullIlu
                }
                Self::AsmIlu0Overlap { .. }
                | Self::RasIlu0Overlap { .. }
                | Self::RasIlukOverlap { .. } => PcDispatchBranch::OverlapIlu,
                Self::IlutLocal => PcDispatchBranch::IlutLocal,
            }
        }

        fn explicit_alias_of(self) -> Option<Self> {
            match self {
                Self::MpiBlockJacobiIlu0Local => Some(Self::Ilu0Local),
                Self::None
                | Self::JacobiWeak
                | Self::Ilu0Local
                | Self::IlutLocal
                | Self::LocalIluk { .. }
                | Self::AsmIlu0Overlap { .. }
                | Self::RasIlu0Overlap { .. }
                | Self::RasIlukOverlap { .. }
                | Self::ReplicatedFullIlu0
                | Self::ReplicatedFullIluk { .. } => None,
            }
        }

        fn semantic_experiment_key(self, mpi_mode: bool) -> String {
            match self {
                Self::MpiBlockJacobiIlu0Local | Self::Ilu0Local if mpi_mode => {
                    "block-jacobi-ilu0-overlap0".to_string()
                }
                Self::MpiBlockJacobiIlu0Local => "mpi-block-jacobi-ilu0-local".to_string(),
                Self::Ilu0Local => "local-ilu0".to_string(),
                Self::None => "none".to_string(),
                Self::JacobiWeak => "jacobi-weak".to_string(),
                Self::IlutLocal => "local-ilut-real-projection-fallback".to_string(),
                Self::LocalIluk { k } => format!("local-iluk:{k}"),
                Self::AsmIlu0Overlap { overlap } => format!("asm-ilu0-overlap{overlap}"),
                Self::RasIlu0Overlap { overlap } => format!("ras-ilu0-overlap{overlap}"),
                Self::RasIlukOverlap { k, overlap } => format!("ras-iluk{k}-overlap{overlap}"),
                Self::ReplicatedFullIlu0 => "replicated-ilu0".to_string(),
                Self::ReplicatedFullIluk { k } => format!("replicated-iluk:{k}"),
            }
        }
        fn is_ilut(self) -> bool {
            matches!(self, Self::IlutLocal)
        }

        fn requires_flexible_right_preconditioning(self) -> bool {
            // Current demo PCs are fixed linear operators. Future variable/nonlinear
            // modes should return true here so run construction preserves FGMRES.
            false
        }
    }

    fn variant_label(variant: FgmresVariant) -> &'static str {
        match variant {
            FgmresVariant::Classical => "classical",
            FgmresVariant::Pipelined => "pipelined",
        }
    }

    fn plain_gmres_variant_label(variant: PlainGmresVariant) -> &'static str {
        match variant {
            PlainGmresVariant::Classical => "classical",
            PlainGmresVariant::Pipelined => "pipelined",
            PlainGmresVariant::SStep { .. } => "s-step",
        }
    }

    fn dist_policy_reason(diag: &DistributedPlanDiagnostics, comm_size: usize) -> &'static str {
        let halo_volume = (diag.halo_recv_volume + diag.halo_send_volume) as f64;
        let overlap_enabled =
            diag.overlap_mode == kryst::matrix::dist_csr::HaloOverlapMode::Interior;
        let communication_heavy = comm_size > 1
            && (overlap_enabled
                || diag.expected_communication_fraction >= 0.55
                || halo_volume >= 4096.0);

        if communication_heavy {
            "high halo/communication pressure detected"
        } else {
            "low halo pressure and compute-dominant local work"
        }
    }

    fn format_dist_policy_report(
        cfg: &BenchmarkConfig,
        sample_spec: Option<&RunSpec>,
        problem: &Problem,
    ) -> String {
        match cfg.dist_policy {
            DistPolicyMode::Off => format!("DistCSR policy: {}", cfg.dist_policy.label()),
            DistPolicyMode::Auto => match sample_spec.map(|spec| spec.ksp).unwrap_or(cfg.ksp) {
                KspKind::Gmres => {
                    let solver = sample_spec
                        .map(|spec| configured_gmres_solver(spec, cfg))
                        .unwrap_or_else(|| GmresSolver::new(1, cfg.rtol, cfg.maxits));
                    format!(
                        "DistCSR policy: auto, KSP=GMRES, selected restart={}, variant={}, reason={}",
                        solver.restart,
                        plain_gmres_variant_label(solver.variant),
                        dist_policy_reason(&problem.dist_plan_diagnostics, problem.comm.size())
                    )
                }
                KspKind::Fgmres => {
                    let mut solver = sample_spec
                        .map(|spec| configured_fgmres_solver(spec, cfg))
                        .unwrap_or_else(|| FgmresSolver::new(cfg.rtol, cfg.maxits, 1));
                    apply_fgmres_dist_plan_policy(&mut solver, problem);
                    format!(
                        "DistCSR policy: auto, KSP=FGMRES, selected restart={}, variant={}, reason={}",
                        solver.restart,
                        variant_label(solver.variant),
                        dist_policy_reason(&problem.dist_plan_diagnostics, problem.comm.size())
                    )
                }
            },
        }
    }

    fn orthog_label(orthog: OrthogMethod) -> &'static str {
        match orthog {
            OrthogMethod::ClassicalGS => "classical-gs",
            OrthogMethod::ModifiedGS => "modified-gs",
        }
    }

    fn reorth_label(reorth: ReorthPolicy) -> &'static str {
        match reorth {
            ReorthPolicy::Never => "never",
            ReorthPolicy::IfNeeded => "if-needed",
            ReorthPolicy::Always => "always",
        }
    }

    fn residual_check_policy_label(policy: ResidualCheckPolicy) -> &'static str {
        match policy {
            ResidualCheckPolicy::RestartOnly => "restart-only",
            ResidualCheckPolicy::OnConvergence => "on-convergence",
            ResidualCheckPolicy::EveryIteration => "every-iteration",
            ResidualCheckPolicy::Debug => "debug",
        }
    }

    fn pc_side_label(pc_side: PcSide) -> &'static str {
        match pc_side {
            PcSide::Right => "right",
            PcSide::Left => "left",
            PcSide::Symmetric => "symmetric",
        }
    }

    fn normalized_fgmres_side(requested_side: PcSide) -> PcSide {
        match requested_side {
            PcSide::Right => PcSide::Right,
            PcSide::Left | PcSide::Symmetric => PcSide::Right,
        }
    }

    #[derive(Clone, Debug)]
    struct ResidualHistoryEntry {
        iter: usize,
        recurrence_residual: R,
        true_residual: Option<R>,
        checkpoint: bool,
    }

    #[derive(Default)]
    struct RunResidualHistory {
        entries: Vec<ResidualHistoryEntry>,
    }

    #[derive(Clone, Debug)]
    struct GlobalReferenceCheck {
        true_residual: R,
        true_residual_rel: R,
        x_error_rel: Option<R>,
    }

    struct RowScaledOp {
        base: Arc<dyn KLinOp<Scalar = S>>,
        d: Vec<R>,
    }

    impl KLinOp for RowScaledOp {
        type Scalar = S;

        fn dims(&self) -> (usize, usize) {
            self.base.dims()
        }

        fn matvec_s(&self, x: &[S], y: &mut [S], scratch: &mut BridgeScratch) {
            self.base.matvec_s(x, y, scratch);
            for (yi, di) in y.iter_mut().zip(self.d.iter().copied()) {
                *yi *= S::from_real(di);
            }
        }

        fn supports_t_matvec_s(&self) -> bool {
            self.base.supports_t_matvec_s()
        }

        fn t_matvec_s(&self, x: &[S], y: &mut [S], scratch: &mut BridgeScratch) {
            self.base.t_matvec_s(x, y, scratch);
        }
    }

    fn compute_row_scaling(matrix: &SparseCsrMatrix<S>, tiny: f64) -> Vec<R> {
        let row_ptr = matrix.row_ptr();
        let vals = matrix.values();
        (0..matrix.nrows())
            .map(|r| {
                let mut row_inf = 0.0f64;
                for nz in row_ptr[r]..row_ptr[r + 1] {
                    row_inf = row_inf.max(vals[nz].abs());
                }
                if row_inf > tiny { 1.0 / row_inf } else { 1.0 }
            })
            .collect()
    }

    fn scale_csr_rows(matrix: &SparseCsrMatrix<S>, row_scale: &[R]) -> SparseCsrMatrix<S> {
        debug_assert_eq!(matrix.nrows(), row_scale.len());
        let row_ptr = matrix.row_ptr().to_vec();
        let col_idx = matrix.col_idx().to_vec();
        let mut values = matrix.values().to_vec();
        for r in 0..matrix.nrows() {
            let d = S::from_real(row_scale[r]);
            for nz in row_ptr[r]..row_ptr[r + 1] {
                values[nz] *= d;
            }
        }
        SparseCsrMatrix::from_csr(matrix.nrows(), matrix.ncols(), row_ptr, col_idx, values)
    }

    fn validate_replicated_full_apply_layout(
        x_local: &[S],
        y_local: &[S],
        global_n: usize,
        global_row_start: usize,
        local_n: usize,
    ) -> Result<(), KError> {
        if x_local.len() != local_n {
            return Err(KError::InvalidInput(format!(
                "replicated ILU apply expected x_local length {local_n}, got {}",
                x_local.len()
            )));
        }
        if y_local.len() != local_n {
            return Err(KError::InvalidInput(format!(
                "replicated ILU apply expected y_local length {local_n}, got {}",
                y_local.len()
            )));
        }
        if global_row_start > global_n || local_n > global_n - global_row_start {
            return Err(KError::InvalidInput(format!(
                "replicated ILU apply invalid owned segment: global_n={global_n}, \
                 global_row_start={global_row_start}, local_n={local_n}"
            )));
        }
        Ok(())
    }

    fn allgather_owned_segment_in_rank_order(
        comm: &UniverseComm,
        x_local: &[S],
        global_n: usize,
        global_row_start: usize,
        local_n: usize,
        x_global: &mut [S],
    ) -> Result<(), KError> {
        validate_replicated_full_apply_layout(
            x_local,
            &vec![S::zero(); local_n],
            global_n,
            global_row_start,
            local_n,
        )?;
        if x_global.len() != global_n {
            return Err(KError::InvalidInput(format!(
                "replicated ILU apply expected x_global scratch length {global_n}, got {}",
                x_global.len()
            )));
        }

        x_global.fill(S::zero());

        #[cfg(feature = "mpi")]
        {
            use mpi::traits::*;

            if let UniverseComm::Mpi(comm_impl) = comm {
                let rank = comm_impl.rank;
                let size = comm_impl.size;
                let local_len = i32::try_from(local_n).map_err(|_| {
                    KError::InvalidInput(format!(
                        "replicated ILU apply local_n={local_n} exceeds MPI i32 length limit"
                    ))
                })?;
                let mut lengths = vec![0i32; size];
                comm_impl
                    .world
                    .all_gather_into(&local_len, &mut lengths[..]);

                let mut offset = 0usize;
                for (r, &len_i32) in lengths.iter().enumerate() {
                    if len_i32 < 0 {
                        return Err(KError::InvalidInput(format!(
                            "replicated ILU apply received negative segment length from rank {r}"
                        )));
                    }
                    let len = len_i32 as usize;
                    if r == rank {
                        if len != x_local.len() || len != local_n {
                            return Err(KError::InvalidInput(format!(
                                "replicated ILU apply local length mismatch: gathered={len}, \
                                 x_local={}, local_n={local_n}",
                                x_local.len()
                            )));
                        }
                        if offset != global_row_start {
                            return Err(KError::InvalidInput(format!(
                                "replicated ILU apply offset mismatch: gathered prefix={offset}, \
                                 global_row_start={global_row_start}"
                            )));
                        }
                    }
                    offset = offset.checked_add(len).ok_or_else(|| {
                        KError::InvalidInput(
                            "replicated ILU apply segment lengths overflowed usize".into(),
                        )
                    })?;
                }
                if offset != global_n {
                    return Err(KError::InvalidInput(format!(
                        "replicated ILU apply expected global_n={global_n}, gathered owned lengths sum={offset}"
                    )));
                }

                let max_len = lengths.iter().copied().max().unwrap_or(0) as usize;
                let mut send = vec![0.0f64; 2 * max_len];
                for (i, value) in x_local.iter().copied().enumerate() {
                    send[2 * i] = value.real();
                    send[2 * i + 1] = value.imag();
                }
                let mut gathered = vec![0.0f64; 2 * max_len * size];
                if max_len > 0 {
                    comm_impl
                        .world
                        .all_gather_into(&send[..], &mut gathered[..]);
                }

                let mut offset = 0usize;
                for (r, &len_i32) in lengths.iter().enumerate() {
                    let len = len_i32 as usize;
                    let start = 2 * r * max_len;
                    for i in 0..len {
                        x_global[offset + i] =
                            S::from_parts(gathered[start + 2 * i], gathered[start + 2 * i + 1]);
                    }
                    offset += len;
                }
                return Ok(());
            }
        }

        if comm.size() != 1 || global_row_start != 0 || local_n != global_n {
            return Err(KError::InvalidInput(
                "replicated ILU apply requires MPI allgather for distributed owned segments".into(),
            ));
        }
        x_global.copy_from_slice(x_local);
        Ok(())
    }

    fn gather_global_reference_vector(
        comm: &UniverseComm,
        local: &[S],
        global_n: usize,
        global_row_start: usize,
        local_n: usize,
        global: &mut [S],
    ) -> Result<(), KError> {
        if local.len() != local_n {
            return Err(KError::InvalidInput(format!(
                "global reference gather expected local length {local_n}, got {}",
                local.len()
            )));
        }
        if global.len() != global_n {
            return Err(KError::InvalidInput(format!(
                "global reference gather expected global scratch length {global_n}, got {}",
                global.len()
            )));
        }
        if global_row_start > global_n || local_n > global_n - global_row_start {
            return Err(KError::InvalidInput(format!(
                "global reference gather invalid owned segment: global_n={global_n}, \
                 global_row_start={global_row_start}, local_n={local_n}"
            )));
        }

        global.fill(S::zero());

        if comm.size() == 1 {
            if global_row_start != 0 || local_n != global_n {
                return Err(KError::InvalidInput(format!(
                    "single-rank global reference gather expected ownership [0..{global_n}), got [{global_row_start}..{})",
                    global_row_start + local_n
                )));
            }
            global.copy_from_slice(local);
            return Ok(());
        }

        #[cfg(feature = "mpi")]
        {
            use mpi::traits::*;

            if let UniverseComm::Mpi(comm_impl) = comm {
                let rank = comm_impl.rank;
                let size = comm_impl.size;
                let local_len = i32::try_from(local_n).map_err(|_| {
                    KError::InvalidInput(format!(
                        "global reference gather local_n={local_n} exceeds MPI i32 length limit"
                    ))
                })?;
                let mut lengths = vec![0i32; size];
                comm_impl
                    .world
                    .all_gather_into(&local_len, &mut lengths[..]);

                let expected_part = DistCsrOp::partition_rows_balanced(global_n, comm);
                let mut offset = 0usize;
                for (r, &len_i32) in lengths.iter().enumerate() {
                    if len_i32 < 0 {
                        return Err(KError::InvalidInput(format!(
                            "global reference gather received negative segment length from rank {r}"
                        )));
                    }
                    let len = len_i32 as usize;
                    let expected_len = expected_part[r + 1] - expected_part[r];
                    if len != expected_len {
                        return Err(KError::InvalidInput(format!(
                            "global reference gather rank {r} length mismatch: gathered={len}, expected balanced={expected_len}"
                        )));
                    }
                    if r == rank {
                        if len != local.len() || len != local_n {
                            return Err(KError::InvalidInput(format!(
                                "global reference gather local length mismatch: gathered={len}, \
                                 local={}, local_n={local_n}",
                                local.len()
                            )));
                        }
                        if offset != global_row_start || expected_part[r] != global_row_start {
                            return Err(KError::InvalidInput(format!(
                                "global reference gather offset mismatch: gathered prefix={offset}, \
                                 expected balanced={}, global_row_start={global_row_start}",
                                expected_part[r]
                            )));
                        }
                    }
                    offset = offset.checked_add(len).ok_or_else(|| {
                        KError::InvalidInput(
                            "global reference gather segment lengths overflowed usize".into(),
                        )
                    })?;
                }
                if offset != global_n {
                    return Err(KError::InvalidInput(format!(
                        "global reference gather expected global_n={global_n}, gathered owned lengths sum={offset}"
                    )));
                }

                let max_len = lengths.iter().copied().max().unwrap_or(0) as usize;
                let mut send = vec![0.0f64; 2 * max_len];
                for (i, value) in local.iter().copied().enumerate() {
                    send[2 * i] = value.real();
                    send[2 * i + 1] = value.imag();
                }
                let mut gathered = vec![0.0f64; 2 * max_len * size];
                if max_len > 0 {
                    comm_impl
                        .world
                        .all_gather_into(&send[..], &mut gathered[..]);
                }

                let mut offset = 0usize;
                for (r, &len_i32) in lengths.iter().enumerate() {
                    let len = len_i32 as usize;
                    let start = 2 * r * max_len;
                    for i in 0..len {
                        global[offset + i] =
                            S::from_parts(gathered[start + 2 * i], gathered[start + 2 * i + 1]);
                    }
                    offset += len;
                }
                return Ok(());
            }
        }

        Err(KError::InvalidInput(
            "global reference gather requires MPI for multi-rank owned segments".into(),
        ))
    }

    fn global_reference_residual(
        problem: &Problem,
        x_local: &[S],
        b_local: &[S],
    ) -> Result<GlobalReferenceCheck, KError> {
        if problem.global_csr.ncols() != problem.global_n {
            return Err(KError::InvalidInput(format!(
                "global reference residual requires square row-owned solution layout, got matrix {}x{} with global_n={}",
                problem.global_csr.nrows(),
                problem.global_csr.ncols(),
                problem.global_n
            )));
        }

        let mut x_global = vec![S::zero(); problem.global_n];
        let mut b_global = vec![S::zero(); problem.global_n];
        gather_global_reference_vector(
            &problem.comm,
            x_local,
            problem.global_n,
            problem.global_row_start,
            problem.local_n,
            &mut x_global,
        )?;
        gather_global_reference_vector(
            &problem.comm,
            b_local,
            problem.global_n,
            problem.global_row_start,
            problem.local_n,
            &mut b_global,
        )?;

        let mut ax_global = vec![S::zero(); problem.global_csr.nrows()];
        problem.global_csr.spmv(&x_global, &mut ax_global);
        let true_global = b_global
            .iter()
            .zip(ax_global.iter())
            .map(|(bi, ai)| (*bi - *ai).abs2())
            .sum::<f64>()
            .sqrt();
        let b_norm = b_global.iter().map(|bi| bi.abs2()).sum::<f64>().sqrt();
        let true_global_rel = true_global / b_norm.max(f64::MIN_POSITIVE);
        let x_error_rel = if problem.rhs_source == RhsSource::GeneratedAOnes
            && problem.solution_reference.valid_for_x_error()
        {
            let err = x_global
                .iter()
                .map(|xi| (*xi - S::one()).abs2())
                .sum::<f64>()
                .sqrt();
            let one_norm = x_global.iter().map(|_| S::one().abs2()).sum::<f64>().sqrt();
            Some(err / one_norm.max(f64::MIN_POSITIVE))
        } else {
            None
        };

        Ok(GlobalReferenceCheck {
            true_residual: true_global,
            true_residual_rel: true_global_rel,
            x_error_rel,
        })
    }

    fn apply_replicated_full_ilu_owned_segment(
        ilu: &IluCsr,
        comm: &UniverseComm,
        side: PcSide,
        x_local: &[S],
        y_local: &mut [S],
        scratch: &mut BridgeScratch,
        global_n: usize,
        global_row_start: usize,
        local_n: usize,
        scratch_in: &std::sync::Mutex<Vec<S>>,
        scratch_out: &std::sync::Mutex<Vec<S>>,
    ) -> Result<(), KError> {
        validate_replicated_full_apply_layout(
            x_local,
            y_local,
            global_n,
            global_row_start,
            local_n,
        )?;

        let mut in_full = scratch_in.lock().map_err(|_| {
            KError::InvalidInput("failed to lock replicated ILU input scratch".into())
        })?;
        let mut out_full = scratch_out.lock().map_err(|_| {
            KError::InvalidInput("failed to lock replicated ILU output scratch".into())
        })?;
        if out_full.len() != global_n {
            return Err(KError::InvalidInput(format!(
                "replicated ILU apply expected output scratch length {global_n}, got {}",
                out_full.len()
            )));
        }

        allgather_owned_segment_in_rank_order(
            comm,
            x_local,
            global_n,
            global_row_start,
            local_n,
            &mut in_full,
        )?;
        out_full.fill(S::zero());
        ilu.apply_s(side, &in_full, &mut out_full, scratch)?;
        y_local.copy_from_slice(&out_full[global_row_start..global_row_start + local_n]);
        Ok(())
    }

    enum PcHandle {
        Jacobi(Jacobi),
        Ilu0(IluCsr),
        MpiBlockJacobiIlu0(IluCsr),
        OverlapIlu(OverlapIluPc),
        ReplicatedFull {
            ilu: IluCsr,
            comm: UniverseComm,
            global_n: usize,
            global_row_start: usize,
            local_n: usize,
            scratch_in: std::sync::Mutex<Vec<S>>,
            scratch_out: std::sync::Mutex<Vec<S>>,
        },
    }

    impl PcHandle {
        fn as_kpc_mut(&mut self) -> &mut dyn KPreconditioner<Scalar = S> {
            match self {
                Self::Jacobi(pc) => pc,
                Self::Ilu0(pc) => pc,
                Self::MpiBlockJacobiIlu0(pc) => pc,
                Self::OverlapIlu(pc) => pc,
                Self::ReplicatedFull { .. } => self,
            }
        }
    }
    impl KPreconditioner for PcHandle {
        type Scalar = S;
        fn dims(&self) -> (usize, usize) {
            match self {
                Self::Jacobi(pc) => KPreconditioner::dims(pc),
                Self::Ilu0(pc) => KPreconditioner::dims(pc),
                Self::MpiBlockJacobiIlu0(pc) => KPreconditioner::dims(pc),
                Self::OverlapIlu(pc) => KPreconditioner::dims(pc),
                Self::ReplicatedFull { global_n, .. } => (*global_n, *global_n),
            }
        }
        fn apply_s(
            &self,
            side: PcSide,
            x: &[S],
            y: &mut [S],
            scratch: &mut BridgeScratch,
        ) -> Result<(), KError> {
            match self {
                Self::Jacobi(pc) => pc.apply_s(side, x, y, scratch),
                Self::Ilu0(pc) => pc.apply_s(side, x, y, scratch),
                Self::MpiBlockJacobiIlu0(pc) => pc.apply_s(side, x, y, scratch),
                Self::OverlapIlu(pc) => pc.apply_s(side, x, y, scratch),
                Self::ReplicatedFull {
                    ilu,
                    comm,
                    global_n,
                    global_row_start,
                    local_n,
                    scratch_in,
                    scratch_out,
                } => apply_replicated_full_ilu_owned_segment(
                    ilu,
                    comm,
                    side,
                    x,
                    y,
                    scratch,
                    *global_n,
                    *global_row_start,
                    *local_n,
                    scratch_in,
                    scratch_out,
                ),
            }
        }
    }

    fn setup_preconditioner_for_run_once(
        problem: &Problem,
        spec: &RunSpec,
        bench_cfg: &BenchmarkConfig,
        pc_csr: &Arc<SparseCsrMatrix<S>>,
        global_pc_csr: &Arc<SparseCsrMatrix<S>>,
        csr_pc_diag: &CsrForPcDiagnostics,
    ) -> Result<Option<PcHandle>, KError> {
        let mut pc: Option<PcHandle> = None;
        match spec.pc.dispatch_branch() {
            PcDispatchBranch::None => {}
            PcDispatchBranch::JacobiWeak => {
                let mut jacobi = Jacobi::new();
                jacobi.setup(pc_csr.as_ref())?;
                pc = Some(PcHandle::Jacobi(jacobi));
            }
            PcDispatchBranch::Ilu0Local => {
                validate_local_ilu_owned_block(pc_csr.as_ref())?;
                let mut cfg = IluCsrConfig::default();
                cfg.kind = match spec.pc {
                    PcKind::Ilu0Local | PcKind::MpiBlockJacobiIlu0Local => IluKind::Ilu0,
                    PcKind::LocalIluk { k } => IluKind::Iluk { k },
                    PcKind::AsmIlu0Overlap { .. }
                    | PcKind::RasIlu0Overlap { .. }
                    | PcKind::RasIlukOverlap { .. }
                    | PcKind::IlutLocal
                    | PcKind::ReplicatedFullIlu0
                    | PcKind::ReplicatedFullIluk { .. }
                    | PcKind::None
                    | PcKind::JacobiWeak => unreachable!(),
                };
                cfg.reordering = bench_cfg.ilu_reordering.clone();
                let mut ilu = IluCsr::new_with_config(cfg);
                ilu.setup(pc_csr.as_ref())?;
                if bench_cfg.run_mode == RunMode::Correctness {
                    let pivot_perturbation_count = None::<usize>;
                    let diag_global = reduce_csr_for_pc_diagnostics(
                        &problem.comm,
                        &csr_pc_diag,
                        pivot_perturbation_count,
                    );
                    if problem.comm.rank() == 0 {
                        print_csr_for_pc_diagnostics(&spec.method_label(), &diag_global);
                    }
                }
                pc = Some(match spec.pc {
                    PcKind::Ilu0Local | PcKind::LocalIluk { .. } => PcHandle::Ilu0(ilu),
                    PcKind::MpiBlockJacobiIlu0Local => PcHandle::MpiBlockJacobiIlu0(ilu),
                    PcKind::AsmIlu0Overlap { .. }
                    | PcKind::RasIlu0Overlap { .. }
                    | PcKind::RasIlukOverlap { .. }
                    | PcKind::IlutLocal
                    | PcKind::ReplicatedFullIlu0
                    | PcKind::ReplicatedFullIluk { .. }
                    | PcKind::None
                    | PcKind::JacobiWeak => unreachable!(),
                });
            }
            PcDispatchBranch::ReplicatedFullIlu => {
                let mut cfg = IluCsrConfig::default();
                cfg.kind = match spec.pc {
                    PcKind::ReplicatedFullIlu0 => IluKind::Ilu0,
                    PcKind::ReplicatedFullIluk { k } => IluKind::Iluk { k },
                    _ => unreachable!(),
                };
                cfg.reordering = bench_cfg.ilu_reordering.clone();
                let mut ilu = IluCsr::new_with_config(cfg);
                ilu.setup(global_pc_csr.as_ref())?;
                if bench_cfg.run_mode == RunMode::Correctness {
                    let pivot_perturbation_count = None::<usize>;
                    let diag_global = reduce_csr_for_pc_diagnostics(
                        &problem.comm,
                        &csr_pc_diag,
                        pivot_perturbation_count,
                    );
                    if problem.comm.rank() == 0 {
                        print_csr_for_pc_diagnostics(&spec.method_label(), &diag_global);
                    }
                }
                pc = Some(PcHandle::ReplicatedFull {
                    ilu,
                    comm: problem.comm.clone(),
                    global_n: problem.global_n,
                    global_row_start: problem.global_row_start,
                    local_n: problem.local_n,
                    scratch_in: std::sync::Mutex::new(vec![S::zero(); problem.global_n]),
                    scratch_out: std::sync::Mutex::new(vec![S::zero(); problem.global_n]),
                });
            }
            PcDispatchBranch::OverlapIlu => {
                let (kind, overlap, restriction) = match spec.pc {
                    PcKind::AsmIlu0Overlap { overlap } => {
                        (IluKind::Ilu0, overlap, OverlapRestriction::Asm)
                    }
                    PcKind::RasIlu0Overlap { overlap } => {
                        (IluKind::Ilu0, overlap, OverlapRestriction::Ras)
                    }
                    PcKind::RasIlukOverlap { k, overlap } => {
                        (IluKind::Iluk { k }, overlap, OverlapRestriction::Ras)
                    }
                    _ => unreachable!(),
                };
                let mut cfg = IluCsrConfig::default();
                cfg.reordering = bench_cfg.ilu_reordering.clone();
                let pc_overlap = OverlapIluPc::setup_from_global_csr(
                    global_pc_csr.as_ref(),
                    problem.comm.clone(),
                    problem.global_row_start,
                    problem.global_row_start + problem.local_n,
                    overlap,
                    kind,
                    cfg,
                    restriction,
                )?;
                pc = Some(PcHandle::OverlapIlu(pc_overlap));
            }
            PcDispatchBranch::IlutLocal => {
                validate_local_ilu_owned_block(pc_csr.as_ref())?;
                let mut cfg = IluCsrConfig::default();
                cfg.kind = IluKind::Ilut {
                    params: Default::default(),
                };
                cfg.reordering = bench_cfg.ilu_reordering.clone();
                let mut ilu = IluCsr::new_with_config(cfg);
                ilu.setup(pc_csr.as_ref())?;
                if bench_cfg.run_mode == RunMode::Correctness {
                    let pivot_perturbation_count = None::<usize>;
                    let diag_global = reduce_csr_for_pc_diagnostics(
                        &problem.comm,
                        &csr_pc_diag,
                        pivot_perturbation_count,
                    );
                    if problem.comm.rank() == 0 {
                        print_csr_for_pc_diagnostics(&spec.method_label(), &diag_global);
                    }
                }
                pc = Some(PcHandle::Ilu0(ilu));
            }
        }
        Ok(pc)
    }

    fn run_once(
        problem: &Problem,
        spec: &RunSpec,
        bench_cfg: &BenchmarkConfig,
    ) -> Result<ResultRow, KError> {
        let b_unscaled = &problem.rhs;
        let (row_scaling, op_scaled): (Option<Vec<R>>, Arc<dyn KLinOp<Scalar = S>>) =
            if bench_cfg.row_scale {
                let d = compute_row_scaling(problem.csr_for_pc.as_ref(), bench_cfg.row_scale_tiny);
                (
                    Some(d.clone()),
                    Arc::new(RowScaledOp {
                        base: problem.op.clone(),
                        d,
                    }),
                )
            } else {
                (None, problem.op.clone())
            };
        let pc_csr: Arc<SparseCsrMatrix<S>> = if bench_cfg.row_scale {
            let row_scale = row_scaling.as_ref().ok_or_else(|| {
                KError::InvalidInput("row scaling requested but factors were not computed".into())
            })?;
            Arc::new(scale_csr_rows(problem.csr_for_pc.as_ref(), row_scale))
        } else {
            problem.csr_for_pc.clone()
        };
        let global_pc_csr: Arc<SparseCsrMatrix<S>> = if bench_cfg.row_scale {
            let row_scale = row_scaling.as_ref().ok_or_else(|| {
                KError::InvalidInput("row scaling requested but factors were not computed".into())
            })?;
            Arc::new(scale_csr_rows(problem.global_csr.as_ref(), row_scale))
        } else {
            problem.global_csr.clone()
        };
        let b_scaled: Vec<S> = if let Some(d) = &row_scaling {
            b_unscaled
                .iter()
                .zip(d.iter())
                .map(|(bi, di)| *bi * S::from_real(*di))
                .collect()
        } else {
            b_unscaled.clone()
        };
        let b = &b_scaled;
        let effective_pc_side = normalized_fgmres_side(spec.pc_side);
        let csr_pc_diag = csr_for_pc_diagnostics(
            pc_csr.as_ref(),
            problem.local_rows_nnz,
            problem.zero_global_rows_local,
        );
        let setup_start = Instant::now();
        let mut pc = setup_preconditioner_for_run_once(
            problem,
            spec,
            bench_cfg,
            &pc_csr,
            &global_pc_csr,
            &csr_pc_diag,
        )?;
        let setup_secs = setup_start.elapsed().as_secs_f64();
        for _ in 0..bench_cfg.warmup_runs {
            let mut x = vec![S::zero(); problem.local_n];
            let mut workspace = Workspace::new(problem.local_n);
            let _ = solve_with_selected_ksp(
                spec,
                bench_cfg,
                problem,
                op_scaled.as_ref(),
                pc.as_mut().map(PcHandle::as_kpc_mut),
                b,
                &mut x,
                effective_pc_side,
                None,
                Some(&mut workspace),
            )?;
        }

        let mut solve_times = Vec::with_capacity(bench_cfg.measured_runs);
        let mut x_last = vec![S::zero(); problem.local_n];
        let mut final_stats = None;
        let mut residual_history_last = RunResidualHistory::default();
        for _ in 0..bench_cfg.measured_runs {
            let mut x = vec![S::zero(); problem.local_n];
            problem.comm.barrier();
            let start = Instant::now();
            let mut workspace = Workspace::new(problem.local_n);
            let mut run_history = RunResidualHistory::default();
            let mut monitors: Vec<Box<MonitorCallback<R>>> = Vec::new();
            if bench_cfg.residual_history {
                monitors.push(Box::new(|it, res, _| {
                    let _ = (it, res);
                    MonitorAction::Continue
                }));
            }
            let history_ref = std::sync::Arc::new(std::sync::Mutex::new(Vec::<(usize, R)>::new()));
            if bench_cfg.residual_history {
                let history_ref_c = history_ref.clone();
                monitors.clear();
                monitors.push(Box::new(move |it, res, _| {
                    if let Ok(mut h) = history_ref_c.lock() {
                        h.push((it, res));
                    }
                    MonitorAction::Continue
                }));
            }
            let stats = solve_with_selected_ksp(
                spec,
                bench_cfg,
                problem,
                op_scaled.as_ref(),
                pc.as_mut().map(PcHandle::as_kpc_mut),
                b,
                &mut x,
                effective_pc_side,
                if monitors.is_empty() {
                    None
                } else {
                    Some(&monitors)
                },
                Some(&mut workspace),
            )?;
            if bench_cfg.residual_history {
                if let Ok(h) = history_ref.lock() {
                    run_history.entries = h
                        .iter()
                        .map(|(it, res)| ResidualHistoryEntry {
                            iter: *it,
                            recurrence_residual: *res,
                            true_residual: None,
                            checkpoint: false,
                        })
                        .collect();
                }
                residual_history_last = run_history;
            }
            problem.comm.barrier();
            let solve_secs = start.elapsed().as_secs_f64();
            solve_times.push(solve_secs);
            x_last = x;
            final_stats = Some(stats);
        }
        let stats = final_stats
            .ok_or_else(|| KError::InvalidInput("no measured solve run executed".into()))?;
        // Row scaling changes only equations (D_r A x = D_r b), so x is unchanged.
        // Keep an explicit "map-back" step so optional future column scaling can hook here.
        let x_unscaled = x_last;
        let min_solve_secs = solve_times.iter().copied().fold(f64::INFINITY, f64::min);
        let median_solve_secs = median(&mut solve_times);

        let reductions = stats.counters.num_global_reductions;
        let overlapped_reduction_waits = stats.counters.overlap_global_reductions;
        let model_predicted_reductions = stats
            .reduction_model
            .as_ref()
            .map(|model| model.estimate_total(stats.iterations));
        let (explicit_true_residual, explicit_true_residual_rel) =
            if bench_cfg.run_mode == RunMode::Correctness {
                let mut ax = vec![S::zero(); b.len()];
                let mut scratch = BridgeScratch::default();
                problem.op.matvec_s(&x_unscaled, &mut ax, &mut scratch);
                for (ri, bi) in ax.iter_mut().zip(b_unscaled.iter().copied()) {
                    *ri = bi - *ri;
                }
                let r2_local = ax.iter().map(|v| v.abs2()).sum::<f64>();
                let true_res = problem.comm.all_reduce_f64(r2_local).sqrt();
                let rhs_norm2_local = b_unscaled.iter().map(|v| v.abs2()).sum::<f64>();
                let rhs_norm = problem.comm.all_reduce_f64(rhs_norm2_local).sqrt();
                let rel_true = true_res / rhs_norm.max(f64::MIN_POSITIVE);
                (Some(true_res), Some(rel_true))
            } else {
                (None, None)
            };
        if bench_cfg.residual_history && bench_cfg.measured_runs > 0 && problem.comm.rank() == 0 {
            let mut checkpoint_count = 0usize;
            // mark restart boundaries based on effective restart interval and policy
            let restart = stats.effective_restart.unwrap_or(spec.restart).max(1);
            for e in &mut residual_history_last.entries {
                if e.iter > 0 && e.iter % restart == 0 {
                    e.checkpoint = true;
                    if matches!(
                        spec.residual_check_policy,
                        ResidualCheckPolicy::RestartOnly
                            | ResidualCheckPolicy::OnConvergence
                            | ResidualCheckPolicy::EveryIteration
                            | ResidualCheckPolicy::Debug
                    ) {
                        e.true_residual = explicit_true_residual;
                    }
                    checkpoint_count += 1;
                }
            }
            println!(
                "[history][rank0] {}: {} points, {} restart checkpoints",
                format!(
                    "{} [row-scale={}]",
                    spec.method_label(),
                    if bench_cfg.row_scale { "on" } else { "off" }
                ),
                residual_history_last.entries.len(),
                checkpoint_count
            );
            if let Some(path) = &bench_cfg.residual_history_file {
                dump_residual_history(path, &residual_history_last)?;
                println!("[history][rank0] wrote {}", path.display());
            }
        }

        let x_error_rel = if bench_cfg.run_mode == RunMode::Correctness
            && problem.rhs_source == RhsSource::GeneratedAOnes
            && problem.solution_reference.valid_for_x_error()
        {
            let err2_local = x_unscaled
                .iter()
                .map(|xi| (*xi - S::one()).abs2())
                .sum::<f64>();
            let one2_local = x_unscaled.iter().map(|_| S::one().abs2()).sum::<f64>();
            let err = problem.comm.all_reduce_f64(err2_local).sqrt();
            let one_norm = problem.comm.all_reduce_f64(one2_local).sqrt();
            Some(err / one_norm.max(f64::MIN_POSITIVE))
        } else {
            None
        };
        let global_reference_check = if bench_cfg.run_mode == RunMode::Correctness {
            Some(global_reference_residual(problem, &x_unscaled, b_unscaled)?)
        } else {
            None
        };
        let verdict_tol = correctness_verdict_tolerance(bench_cfg);
        let dist_ok = explicit_true_residual_rel.map(|rel| rel <= verdict_tol);
        let global_ok = global_reference_check
            .as_ref()
            .map(|check| check.true_residual_rel <= verdict_tol);
        let x_ok = if bench_cfg.run_mode == RunMode::Correctness
            && problem.rhs_source == RhsSource::GeneratedAOnes
            && problem.solution_reference.valid_for_x_error()
        {
            global_reference_check
                .as_ref()
                .and_then(|check| check.x_error_rel)
                .map(|rel| rel <= verdict_tol)
        } else {
            None
        };
        let dof_per_sec = if matches!(
            problem.backend,
            CsrBackend::Serial | CsrBackend::Distributed
        ) && median_solve_secs > 0.0
        {
            Some(problem.global_n as f64 / median_solve_secs)
        } else {
            None
        };

        Ok(ResultRow {
            operator_storage: operator_storage_label(problem),
            execution_backend: execution_backend_label(problem),
            pc_domain: pc_domain_label(spec.pc, problem),
            pc_apply: pc_apply_label(spec.pc, problem),
            method: format!(
                "{} [row-scale={}]",
                spec.method_label(),
                if bench_cfg.row_scale { "on" } else { "off" }
            ),
            requested_policy: spec.requested_policy_label(),
            effective_policy: format!(
                "ksp={}, variant={}, restart={}, residual-check={}",
                spec.ksp.label(),
                stats
                    .effective_variant
                    .as_deref()
                    .unwrap_or(variant_label(spec.variant)),
                stats.effective_restart.unwrap_or(spec.restart),
                stats
                    .effective_residual_check_policy
                    .as_deref()
                    .unwrap_or(residual_check_policy_label(spec.residual_check_policy))
            ) + if bench_cfg.run_mode == RunMode::Correctness
                && bench_cfg.mark_replicated_check
            {
                " [replicated-check=enabled]"
            } else {
                ""
            },
            setup_secs,
            median_solve_secs,
            min_solve_secs,
            iterations: stats.iterations,
            reductions,
            overlapped_reduction_waits,
            model_predicted_reductions,
            restart_count: stats.fgmres_counters.as_ref().map(|c| c.restart_count),
            inner_iterations_last_cycle: stats
                .fgmres_counters
                .as_ref()
                .map(|c| c.inner_iterations_last_cycle),
            pipeline_fallbacks: stats.fgmres_counters.as_ref().map(|c| c.pipeline_fallbacks),
            reported_residual: stats.final_residual,
            explicit_true_residual,
            explicit_true_residual_rel,
            x_error_rel,
            global_true_residual: global_reference_check
                .as_ref()
                .map(|check| check.true_residual),
            global_true_residual_rel: global_reference_check
                .as_ref()
                .map(|check| check.true_residual_rel),
            global_x_error_rel: global_reference_check.and_then(|check| check.x_error_rel),
            dist_ok,
            global_ok,
            x_ok,
            reason: stats.reason,
            dof_per_sec,
        })
    }

    fn validate_local_ilu_owned_block(matrix: &SparseCsrMatrix<S>) -> Result<(), KError> {
        if matrix.nrows() != matrix.ncols() {
            return Err(KError::InvalidInput(
                "local ILU requires square owned block".into(),
            ));
        }
        Ok(())
    }

    fn csr_for_pc_diagnostics(
        matrix: &SparseCsrMatrix<S>,
        nnz_local_rows: usize,
        zero_global_rows: usize,
    ) -> CsrForPcDiagnostics {
        let mut diag_min_abs = f64::INFINITY;
        let mut diag_max_abs = 0.0f64;
        let mut zero_local_rows = 0usize;
        let mut structural_diag_missing = 0usize;
        let mut diag_tiny_or_missing_count = 0usize;
        let tiny_diag_threshold = 1e-14f64;
        let row_ptr = matrix.row_ptr();
        let col_idx = matrix.col_idx();
        let values = matrix.values();
        for r in 0..matrix.nrows() {
            if row_ptr[r] == row_ptr[r + 1] {
                zero_local_rows += 1;
            }
            let mut diag = None;
            for nz in row_ptr[r]..row_ptr[r + 1] {
                if col_idx[nz] == r {
                    diag = Some(values[nz].abs());
                    break;
                }
            }
            match diag {
                Some(v) => {
                    diag_min_abs = diag_min_abs.min(v);
                    diag_max_abs = diag_max_abs.max(v);
                    if v <= tiny_diag_threshold {
                        diag_tiny_or_missing_count += 1;
                    }
                }
                None => {
                    structural_diag_missing += 1;
                    diag_tiny_or_missing_count += 1;
                }
            }
        }
        if !diag_min_abs.is_finite() {
            diag_min_abs = 0.0;
        }
        let nnz_local_block = values.len();
        let nnz_ratio = nnz_local_block as f64 / (nnz_local_rows.max(1) as f64);
        CsrForPcDiagnostics {
            nnz_local_block,
            nnz_local_rows,
            nnz_ratio,
            zero_local_rows,
            zero_global_rows,
            structural_diag_missing,
            diag_min_abs,
            diag_max_abs,
            diag_tiny_or_missing_count,
        }
    }

    fn rank_spread_from_locals(locals: &[f64], rank0_local: f64) -> RankSpread {
        let global_sum = locals.iter().sum::<f64>();
        let rank_min = locals.iter().copied().fold(f64::INFINITY, f64::min);
        let rank_max = locals.iter().copied().fold(f64::NEG_INFINITY, f64::max);
        RankSpread {
            global_sum,
            rank_min: if rank_min.is_finite() { rank_min } else { 0.0 },
            rank_max: if rank_max.is_finite() { rank_max } else { 0.0 },
            rank0_local,
        }
    }

    fn reduce_rank_spread<C: Comm<Vec = Vec<f64>>>(comm: &C, local: f64) -> RankSpread {
        if comm.size() <= 1 {
            return RankSpread {
                global_sum: local,
                rank_min: local,
                rank_max: local,
                rank0_local: local,
            };
        }
        let global_sum = comm.all_reduce_f64(local);
        let mut gathered = Vec::new();
        comm.gather(&[local], &mut gathered, 0);
        let (rank_min, rank_max) = if comm.rank() == 0 && !gathered.is_empty() {
            (
                gathered.iter().copied().fold(f64::INFINITY, f64::min),
                gathered.iter().copied().fold(f64::NEG_INFINITY, f64::max),
            )
        } else {
            (0.0, 0.0)
        };
        let rank0_local = if comm.rank() == 0 { local } else { 0.0 };
        RankSpread {
            global_sum,
            rank_min,
            rank_max,
            rank0_local,
        }
    }

    fn reduce_csr_for_pc_diagnostics<C: Comm<Vec = Vec<f64>>>(
        comm: &C,
        local: &CsrForPcDiagnostics,
        pivot_perturbation_count: Option<usize>,
    ) -> CsrForPcDiagnosticsGlobal {
        let zero_local_rows_rank_count =
            reduce_rank_spread(comm, if local.zero_local_rows > 0 { 1.0 } else { 0.0 }).global_sum;
        CsrForPcDiagnosticsGlobal {
            nnz_local_block: reduce_rank_spread(comm, local.nnz_local_block as f64),
            nnz_local_rows: reduce_rank_spread(comm, local.nnz_local_rows as f64),
            nnz_ratio: reduce_rank_spread(comm, local.nnz_ratio),
            zero_local_rows_rank_count,
            zero_global_rows: reduce_rank_spread(comm, local.zero_global_rows as f64),
            structural_diag_missing: reduce_rank_spread(comm, local.structural_diag_missing as f64),
            diag_min_abs: reduce_rank_spread(comm, local.diag_min_abs),
            diag_max_abs: reduce_rank_spread(comm, local.diag_max_abs),
            diag_tiny_or_missing_count: reduce_rank_spread(
                comm,
                local.diag_tiny_or_missing_count as f64,
            ),
            pivot_perturbation_count: pivot_perturbation_count
                .map(|v| reduce_rank_spread(comm, v as f64)),
        }
    }

    fn format_csr_for_pc_diagnostics(
        method_label: &str,
        diag: &CsrForPcDiagnosticsGlobal,
    ) -> String {
        let pivot_text = diag
            .pivot_perturbation_count
            .as_ref()
            .map(|v| {
                format!(
                    "global_sum={:.0}, rank_min={:.0}, rank_max={:.0}, rank0_local={:.0}",
                    v.global_sum, v.rank_min, v.rank_max, v.rank0_local
                )
            })
            .unwrap_or_else(|| "N/A".to_string());
        let nnz_ratio_global =
            diag.nnz_local_block.global_sum / diag.nnz_local_rows.global_sum.max(1.0);
        format!(
            "[diag][rank0] {method_label}: nnz(A_ii):global_sum={:.0},rank_min={:.0},rank_max={:.0},rank0_local={:.0}; nnz(local rows):global_sum={:.0},rank_min={:.0},rank_max={:.0},rank0_local={:.0}; nnz_ratio_global={:.3e},rank_local_min={:.3e},rank_local_max={:.3e},rank0_local={:.3e}; zero_local_rows_rank_count={:.0}; zero_global_rows:global_sum={:.0},rank_min={:.0},rank_max={:.0},rank0_local={:.0}; structural_diag_missing:global_sum={:.0},rank_min={:.0},rank_max={:.0},rank0_local={:.0}; |diag|min:global_sum={:.3e},rank_min={:.3e},rank_max={:.3e},rank0_local={:.3e}; |diag|max:global_sum={:.3e},rank_min={:.3e},rank_max={:.3e},rank0_local={:.3e}; tiny/missing diag:global_sum={:.0},rank_min={:.0},rank_max={:.0},rank0_local={:.0}; pivot perturbations={}",
            diag.nnz_local_block.global_sum,
            diag.nnz_local_block.rank_min,
            diag.nnz_local_block.rank_max,
            diag.nnz_local_block.rank0_local,
            diag.nnz_local_rows.global_sum,
            diag.nnz_local_rows.rank_min,
            diag.nnz_local_rows.rank_max,
            diag.nnz_local_rows.rank0_local,
            nnz_ratio_global,
            diag.nnz_ratio.rank_min,
            diag.nnz_ratio.rank_max,
            diag.nnz_ratio.rank0_local,
            diag.zero_local_rows_rank_count,
            diag.zero_global_rows.global_sum,
            diag.zero_global_rows.rank_min,
            diag.zero_global_rows.rank_max,
            diag.zero_global_rows.rank0_local,
            diag.structural_diag_missing.global_sum,
            diag.structural_diag_missing.rank_min,
            diag.structural_diag_missing.rank_max,
            diag.structural_diag_missing.rank0_local,
            diag.diag_min_abs.global_sum,
            diag.diag_min_abs.rank_min,
            diag.diag_min_abs.rank_max,
            diag.diag_min_abs.rank0_local,
            diag.diag_max_abs.global_sum,
            diag.diag_max_abs.rank_min,
            diag.diag_max_abs.rank_max,
            diag.diag_max_abs.rank0_local,
            diag.diag_tiny_or_missing_count.global_sum,
            diag.diag_tiny_or_missing_count.rank_min,
            diag.diag_tiny_or_missing_count.rank_max,
            diag.diag_tiny_or_missing_count.rank0_local,
            pivot_text
        )
    }

    fn print_csr_for_pc_diagnostics(method_label: &str, diag: &CsrForPcDiagnosticsGlobal) {
        println!("{}", format_csr_for_pc_diagnostics(method_label, diag));
    }

    fn configured_fgmres_solver(spec: &RunSpec, bench_cfg: &BenchmarkConfig) -> FgmresSolver {
        let mut solver = FgmresSolver::new(bench_cfg.rtol, bench_cfg.maxits, spec.restart);
        solver.variant = spec.variant;
        solver.residual_check_policy = spec.residual_check_policy;
        solver.orthog = spec.orthog;
        solver.reorth = spec.reorth;
        solver.atol = bench_cfg.atol;
        solver.dtol = 1e6;
        solver.haptol = bench_cfg.fgmres_haptol;
        solver.min_inner_before_fallback = bench_cfg.min_inner_before_fallback.max(1);
        solver.stagnation_policy =
            if bench_cfg.run_mode == RunMode::Correctness || !bench_cfg.allow_stagnation_fallback {
                FgmresStagnationPolicy::Disabled
            } else {
                match spec.variant {
                    FgmresVariant::Classical => FgmresStagnationPolicy::RestartClassicalToo,
                    FgmresVariant::Pipelined => FgmresStagnationPolicy::PipelineFallbackOnly,
                }
            };
        solver
    }

    fn apply_fgmres_dist_plan_policy(solver: &mut FgmresSolver, problem: &Problem) {
        solver.apply_distcsr_policy(&problem.dist_plan_diagnostics, problem.comm.size());
    }

    fn gmres_variant_from_fgmres(variant: FgmresVariant) -> PlainGmresVariant {
        match variant {
            FgmresVariant::Classical => PlainGmresVariant::Classical,
            FgmresVariant::Pipelined => PlainGmresVariant::Pipelined,
        }
    }

    fn gmres_orthog_from_fgmres(orthog: OrthogMethod) -> GmresOrthog {
        match orthog {
            OrthogMethod::ClassicalGS => GmresOrthog::Cgs,
            OrthogMethod::ModifiedGS => GmresOrthog::Mgs,
        }
    }

    fn configured_gmres_solver(spec: &RunSpec, bench_cfg: &BenchmarkConfig) -> GmresSolver {
        let mut solver = GmresSolver::new(spec.restart, bench_cfg.rtol, bench_cfg.maxits);
        solver.variant = gmres_variant_from_fgmres(spec.variant);
        solver.orthog = gmres_orthog_from_fgmres(spec.orthog);
        solver.reorth = spec.reorth;
        solver.conv.atol = bench_cfg.atol;
        solver.conv.dtol = 1e6;
        solver.haptol = bench_cfg.fgmres_haptol;
        solver
    }

    #[allow(clippy::too_many_arguments)]
    fn solve_with_selected_ksp(
        spec: &RunSpec,
        bench_cfg: &BenchmarkConfig,
        problem: &Problem,
        op: &(dyn KLinOp<Scalar = S> + '_),
        mut pc: Option<&mut dyn KPreconditioner<Scalar = S>>,
        b: &[S],
        x: &mut [S],
        pc_side: PcSide,
        monitors: Option<&[Box<MonitorCallback<R>>]>,
        mut work: Option<&mut Workspace>,
    ) -> Result<kryst::utils::convergence::SolveStats<R>, KError> {
        match spec.ksp {
            KspKind::Gmres => {
                let mut solver = configured_gmres_solver(spec, bench_cfg);
                if let Some(w) = work.as_deref_mut() {
                    solver.setup_workspace(w);
                }
                let pc_ref = pc
                    .as_deref_mut()
                    .map(|p| p as &dyn KPreconditioner<Scalar = S>);
                solver.solve(op, pc_ref, b, x, pc_side, &problem.comm, monitors, work)
            }
            KspKind::Fgmres => {
                let mut solver = configured_fgmres_solver(spec, bench_cfg);
                if bench_cfg.dist_policy == DistPolicyMode::Auto {
                    apply_fgmres_dist_plan_policy(&mut solver, problem);
                }
                if let Some(w) = work.as_deref_mut() {
                    solver.setup_workspace(w);
                }
                solver.solve_k(op, pc, b, x, pc_side, &problem.comm, monitors, work)
            }
        }
    }

    fn dump_residual_history(path: &Path, history: &RunResidualHistory) -> Result<(), KError> {
        let mut f = File::create(path)
            .map_err(|e| KError::InvalidInput(format!("failed to create history file: {e}")))?;
        if path.extension().and_then(|e| e.to_str()) == Some("json") {
            writeln!(f, "[").map_err(|e| KError::InvalidInput(format!("write failed: {e}")))?;
            for (i, e) in history.entries.iter().enumerate() {
                writeln!(f, "  {{\"iter\":{},\"recurrence_residual\":{},\"true_residual\":{},\"checkpoint\":{}}}{}", e.iter, e.recurrence_residual, e.true_residual.map(|v| v.to_string()).unwrap_or_else(|| "null".to_string()), e.checkpoint, if i + 1 == history.entries.len() {""} else {","})
                    .map_err(|e| KError::InvalidInput(format!("write failed: {e}")))?;
            }
            writeln!(f, "]").map_err(|e| KError::InvalidInput(format!("write failed: {e}")))?;
        } else {
            writeln!(f, "iter,recurrence_residual,true_residual,checkpoint")
                .map_err(|e| KError::InvalidInput(format!("write failed: {e}")))?;
            for e in &history.entries {
                writeln!(
                    f,
                    "{},{},{},{}",
                    e.iter,
                    e.recurrence_residual,
                    e.true_residual.map(|v| v.to_string()).unwrap_or_default(),
                    e.checkpoint
                )
                .map_err(|e| KError::InvalidInput(format!("write failed: {e}")))?;
            }
        }
        Ok(())
    }

    fn median(samples: &mut [f64]) -> f64 {
        samples.sort_by(f64::total_cmp);
        let n = samples.len();
        if n % 2 == 1 {
            samples[n / 2]
        } else {
            (samples[n / 2 - 1] + samples[n / 2]) * 0.5
        }
    }

    fn render_verdict(flag: Option<bool>) -> &'static str {
        match flag {
            Some(true) => "OK",
            Some(false) => "FAIL",
            None => "N/A",
        }
    }

    fn render_reason(row: &ResultRow, problem: &Problem) -> String {
        if row.reason == ConvergedReason::ConvergedRtol
            && problem.generated_case
            && problem.comm.size() > 1
            && (row.global_ok == Some(false) || row.x_ok == Some(false))
        {
            "ConvergedRtolDistOpButFailedGlobalReference".to_string()
        } else {
            format!("{:?}", row.reason)
        }
    }

    fn render_result_row(row: &ResultRow, mode: RunMode, problem: &Problem) -> String {
        let rst = row
            .restart_count
            .map(|v| v.to_string())
            .unwrap_or_else(|| "N/A".to_string());
        let inn = row
            .inner_iterations_last_cycle
            .map(|v| v.to_string())
            .unwrap_or_else(|| "N/A".to_string());
        let pfb = row
            .pipeline_fallbacks
            .map(|v| v.to_string())
            .unwrap_or_else(|| "N/A".to_string());
        let include_dof_col = matches!(
            problem.backend,
            CsrBackend::Serial | CsrBackend::Distributed
        );
        if mode == RunMode::Scalability {
            let dof = row
                .dof_per_sec
                .map(|v| format!("{v:.2e}"))
                .unwrap_or_else(|| "N/A".to_string());
            return format!(
                "{:<7} {:<8} {:<6} {:<42} {:<36} {:<34} {:>9.3} {:>9.3} {:>7} {:>5} {:>5} {:>5} {:>17} {:>14.2e} {:>12}",
                row.operator_storage,
                row.execution_backend,
                row.pc_domain,
                row.pc_apply,
                row.method,
                row.effective_policy,
                row.median_solve_secs,
                row.min_solve_secs,
                row.iterations,
                rst,
                inn,
                pfb,
                format!(
                    "{}/{}/{}",
                    row.reductions,
                    row.overlapped_reduction_waits,
                    row.model_predicted_reductions
                        .map(|v| v.to_string())
                        .unwrap_or_else(|| "n/a".to_string())
                ),
                row.reported_residual,
                dof
            );
        }
        let explicit_true = row
            .explicit_true_residual
            .map(|v| format!("{v:.2e}"))
            .unwrap_or_else(|| "N/A".to_string());
        let explicit_true_rel = row
            .explicit_true_residual_rel
            .map(|v| format!("{v:.2e}"))
            .unwrap_or_else(|| "N/A".to_string());
        let x_error_rel = row
            .x_error_rel
            .map(|v| format!("{v:.2e}"))
            .unwrap_or_else(|| "N/A".to_string());
        let global_true = row
            .global_true_residual
            .map(|v| format!("{v:.2e}"))
            .unwrap_or_else(|| "N/A".to_string());
        let global_true_rel = row
            .global_true_residual_rel
            .map(|v| format!("{v:.2e}"))
            .unwrap_or_else(|| "N/A".to_string());
        let global_x_error_rel = row
            .global_x_error_rel
            .map(|v| format!("{v:.2e}"))
            .unwrap_or_else(|| "N/A".to_string());
        let dist_ok = render_verdict(row.dist_ok);
        let global_ok = render_verdict(row.global_ok);
        let x_ok = render_verdict(row.x_ok);
        let reason = render_reason(row, problem);
        if include_dof_col {
            let dof = row
                .dof_per_sec
                .map(|v| format!("{v:.2e}"))
                .unwrap_or_else(|| "N/A".to_string());
            format!(
                "{:<7} {:<8} {:<6} {:<42} {:<36} {:<34} {:<34} {:>9.3} {:>9.3} {:>9.3} {:>7} {:>5} {:>5} {:>5} {:>14.2e} {:>14} {:>14} {:>14} {:>14} {:>14} {:>14} {:>8} {:>8} {:>8} {:>26} {:>12}",
                row.operator_storage,
                row.execution_backend,
                row.pc_domain,
                row.pc_apply,
                row.method,
                row.requested_policy,
                row.effective_policy,
                row.setup_secs,
                row.median_solve_secs,
                row.min_solve_secs,
                row.iterations,
                rst,
                inn,
                pfb,
                row.reported_residual,
                explicit_true,
                explicit_true_rel,
                x_error_rel,
                global_true,
                global_true_rel,
                global_x_error_rel,
                dist_ok,
                global_ok,
                x_ok,
                reason,
                dof
            )
        } else {
            format!(
                "{:<7} {:<8} {:<6} {:<42} {:<36} {:<34} {:<34} {:>9.3} {:>9.3} {:>9.3} {:>7} {:>5} {:>5} {:>5} {:>14.2e} {:>14} {:>14} {:>14} {:>14} {:>14} {:>14} {:>8} {:>8} {:>8} {:>26}",
                row.operator_storage,
                row.execution_backend,
                row.pc_domain,
                row.pc_apply,
                row.method,
                row.requested_policy,
                row.effective_policy,
                row.setup_secs,
                row.median_solve_secs,
                row.min_solve_secs,
                row.iterations,
                rst,
                inn,
                pfb,
                row.reported_residual,
                explicit_true,
                explicit_true_rel,
                x_error_rel,
                global_true,
                global_true_rel,
                global_x_error_rel,
                dist_ok,
                global_ok,
                x_ok,
                reason
            )
        }
    }

    fn render_failure_result_row(spec: &RunSpec, mode: RunMode, problem: &Problem) -> String {
        let include_dof_col = matches!(
            problem.backend,
            CsrBackend::Serial | CsrBackend::Distributed
        );
        if mode == RunMode::Scalability {
            return format!(
                "{:<7} {:<8} {:<6} {:<42} {:<36} {:<34} {:>9} {:>9} {:>7} {:>5} {:>5} {:>5} {:>17} {:>14} {:>12}",
                operator_storage_label(problem),
                execution_backend_label(problem),
                pc_domain_label(spec.pc, problem),
                pc_apply_label(spec.pc, problem),
                spec.method_label(),
                "N/A",
                "FAIL",
                "FAIL",
                "N/A",
                "N/A",
                "N/A",
                "N/A",
                "N/A",
                "N/A",
                "N/A"
            );
        }
        if include_dof_col {
            format!(
                "{:<7} {:<8} {:<6} {:<42} {:<36} {:<34} {:<34} {:>9} {:>9} {:>9} {:>7} {:>5} {:>5} {:>5} {:>14} {:>14} {:>14} {:>14} {:>14} {:>14} {:>14} {:>8} {:>8} {:>8} {:>26} {:>12}",
                operator_storage_label(problem),
                execution_backend_label(problem),
                pc_domain_label(spec.pc, problem),
                pc_apply_label(spec.pc, problem),
                spec.method_label(),
                spec.requested_policy_label(),
                "N/A",
                "FAIL",
                "FAIL",
                "FAIL",
                "N/A",
                "N/A",
                "N/A",
                "N/A",
                "N/A",
                "N/A",
                "N/A",
                "N/A",
                "N/A",
                "N/A",
                "N/A",
                "N/A",
                "N/A",
                "N/A",
                "N/A",
                "N/A"
            )
        } else {
            format!(
                "{:<7} {:<8} {:<6} {:<42} {:<36} {:<34} {:<34} {:>9} {:>9} {:>9} {:>7} {:>5} {:>5} {:>5} {:>14} {:>14} {:>14} {:>14} {:>14} {:>14} {:>14} {:>8} {:>8} {:>8} {:>26}",
                operator_storage_label(problem),
                execution_backend_label(problem),
                pc_domain_label(spec.pc, problem),
                pc_apply_label(spec.pc, problem),
                spec.method_label(),
                spec.requested_policy_label(),
                "N/A",
                "FAIL",
                "FAIL",
                "FAIL",
                "N/A",
                "N/A",
                "N/A",
                "N/A",
                "N/A",
                "N/A",
                "N/A",
                "N/A",
                "N/A",
                "N/A",
                "N/A",
                "N/A",
                "N/A",
                "N/A",
                "N/A"
            )
        }
    }

    fn solution_reference_diagnostics(matrix: &SparseCsrMatrix<S>) -> SolutionReferenceDiagnostics {
        let tiny_diag_threshold = 1e-14f64;
        let row_ptr = matrix.row_ptr();
        let col_idx = matrix.col_idx();
        let values = matrix.values();
        let mut zero_row_count = 0usize;
        let mut missing_diagonal_count = 0usize;
        let mut tiny_diagonal_count = 0usize;

        for row in 0..matrix.nrows() {
            let start = row_ptr[row];
            let end = row_ptr[row + 1];
            if start == end {
                zero_row_count += 1;
            }
            let mut diag = None;
            for nz in start..end {
                if col_idx[nz] == row {
                    diag = Some(values[nz].abs());
                    break;
                }
            }
            match diag {
                Some(v) if v <= tiny_diag_threshold => tiny_diagonal_count += 1,
                Some(_) => {}
                None => missing_diagonal_count += 1,
            }
        }

        SolutionReferenceDiagnostics {
            global_rows: matrix.nrows(),
            global_cols: matrix.ncols(),
            zero_row_count,
            missing_diagonal_count,
            tiny_diagonal_count,
        }
    }

    fn build_shifted_poisson_case(
        comm: &UniverseComm,
        nx: usize,
        ny: usize,
        alpha: f64,
        beta: f64,
        convection_x: f64,
        convection_y: f64,
    ) -> Result<Problem, KError> {
        if nx == 0 || ny == 0 {
            return Err(KError::InvalidInput(
                "shifted Poisson grid dimensions must be positive".into(),
            ));
        }
        let n = nx.checked_mul(ny).ok_or_else(|| {
            KError::InvalidInput(format!("shifted Poisson grid {nx}x{ny} overflows usize"))
        })?;
        if n < comm.size() {
            return Err(KError::InvalidInput(format!(
                "shifted Poisson grid {nx}x{ny} has {n} rows, fewer than {} MPI ranks",
                comm.size()
            )));
        }

        let csr_sparse = build_shifted_poisson_csr(nx, ny, alpha, beta, convection_x, convection_y);
        let solution_reference = solution_reference_diagnostics(&csr_sparse);
        let row_part = DistCsrOp::partition_rows_balanced(n, comm);
        let row_start = row_part[comm.rank()];
        let row_end = row_part[comm.rank() + 1];
        let local_csr = slice_csr_rows(&csr_sparse, row_start, row_end);
        let local_pc_block = slice_csr_rows_owned_cols(&csr_sparse, row_start, row_end);

        let op = DistCsrOp::from_local_rows(n, row_start, &local_csr, &row_part, comm.clone())?;
        let dist_plan_diagnostics = op.plan_diagnostics().clone();
        let op_arc: Arc<dyn KLinOp<Scalar = S>> = Arc::new(op);

        let ones = vec![S::one(); n];
        let mut rhs_global = vec![S::zero(); n];
        csr_sparse.spmv(&ones, &mut rhs_global);
        let rhs = rhs_global[row_start..row_end].to_vec();

        let local_n = row_end - row_start;
        let backend = classify_backend(comm.size(), local_n, n);

        Ok(Problem {
            op: op_arc,
            dist_plan_diagnostics,
            rhs,
            rhs_source: RhsSource::GeneratedAOnes,
            solution_reference,
            csr_for_pc: Arc::new(local_pc_block),
            global_csr: Arc::new(csr_sparse),
            local_rows_nnz: local_csr.values().len(),
            zero_global_rows_local: count_zero_rows(&local_csr),
            local_n,
            global_n: n,
            global_row_start: row_start,
            comm: comm.clone(),
            backend,
            backend_descr: format!(
                "Generated shifted Poisson CSR ({}x{} grid, shift={alpha:.3}+{beta:.3}i, convection=({convection_x:.3},{convection_y:.3}))",
                nx, ny
            ),
            generated_case: true,
        })
    }

    fn build_shifted_poisson_csr(
        nx: usize,
        ny: usize,
        alpha: f64,
        beta: f64,
        convection_x: f64,
        convection_y: f64,
    ) -> SparseCsrMatrix<S> {
        let n = nx * ny;
        let shift = S::from_parts(alpha, beta);
        let mut row_ptr = Vec::with_capacity(n + 1);
        let mut col_idx = Vec::with_capacity(n * 5);
        let mut values = Vec::with_capacity(n * 5);
        row_ptr.push(0);

        for iy in 0..ny {
            for ix in 0..nx {
                let row = iy * nx + ix;
                let mut entries = Vec::with_capacity(5);
                if iy > 0 {
                    entries.push((row - nx, S::from_real(-1.0 - 0.5 * convection_y)));
                }
                if ix > 0 {
                    entries.push((row - 1, S::from_real(-1.0 - 0.5 * convection_x)));
                }
                entries.push((row, S::from_real(4.0) + shift));
                if ix + 1 < nx {
                    entries.push((row + 1, S::from_real(-1.0 + 0.5 * convection_x)));
                }
                if iy + 1 < ny {
                    entries.push((row + nx, S::from_real(-1.0 + 0.5 * convection_y)));
                }
                entries.sort_by_key(|&(col, _)| col);
                for (col, val) in entries {
                    col_idx.push(col);
                    values.push(val);
                }
                row_ptr.push(col_idx.len());
            }
        }

        SparseCsrMatrix::from_csr(n, n, row_ptr, col_idx, values)
    }

    fn load_problem_complex(
        mat_path: &Path,
        comm: &UniverseComm,
        _run_mode: RunMode,
        _mark_replicated_check: bool,
    ) -> Result<Problem, KError> {
        let mm = read_matrix_market(mat_path)?;
        let csr_sparse: SparseCsrMatrix<S> = mm.to_csr_matrix_scalar()?;
        let nrows = csr_sparse.nrows();
        let ncols = csr_sparse.ncols();
        let solution_reference = solution_reference_diagnostics(&csr_sparse);

        let row_part = DistCsrOp::partition_rows_balanced(nrows, comm);
        let row_start = row_part[comm.rank()];
        let row_end = row_part[comm.rank() + 1];
        let local_csr = slice_csr_rows(&csr_sparse, row_start, row_end);
        let local_pc_block = slice_csr_rows_owned_cols(&csr_sparse, row_start, row_end);

        let op = DistCsrOp::from_local_rows(nrows, row_start, &local_csr, &row_part, comm.clone())?;
        let dist_plan_diagnostics = op.plan_diagnostics().clone();
        let op_arc: Arc<dyn KLinOp<Scalar = S>> = Arc::new(op);

        let (rhs_global, rhs_source) = match try_load_rhs_s(mat_path, nrows) {
            Some(vec) => (vec, RhsSource::LoadedFromFile),
            None => {
                let ones = vec![S::one(); ncols];
                let mut b = vec![S::zero(); nrows];
                csr_sparse.spmv(&ones, &mut b);
                (b, RhsSource::GeneratedAOnes)
            }
        };
        let rhs = rhs_global[row_start..row_end].to_vec();

        let local_n = row_end - row_start;
        let global_n = nrows;
        let backend = classify_backend(comm.size(), local_n, global_n);

        Ok(Problem {
            op: op_arc,
            dist_plan_diagnostics,
            rhs,
            rhs_source,
            solution_reference,
            csr_for_pc: Arc::new(local_pc_block),
            global_csr: Arc::new(csr_sparse),
            local_rows_nnz: local_csr.values().len(),
            zero_global_rows_local: count_zero_rows(&local_csr),
            local_n,
            global_n,
            global_row_start: row_start,
            comm: comm.clone(),
            backend,
            backend_descr: "Distributed CSR (complex)".to_string(),
            generated_case: false,
        })
    }

    fn count_zero_rows(matrix: &SparseCsrMatrix<S>) -> usize {
        let row_ptr = matrix.row_ptr();
        (0..matrix.nrows())
            .filter(|&r| row_ptr[r] == row_ptr[r + 1])
            .count()
    }

    fn slice_csr_rows(matrix: &SparseCsrMatrix<S>, start: usize, end: usize) -> SparseCsrMatrix<S> {
        let row_ptr = matrix.row_ptr();
        let col_idx = matrix.col_idx();
        let values = matrix.values();
        let start_nnz = row_ptr[start];
        let end_nnz = row_ptr[end];

        let mut local_rp = Vec::with_capacity(end - start + 1);
        for r in start..=end {
            local_rp.push(row_ptr[r] - start_nnz);
        }
        let local_ci = col_idx[start_nnz..end_nnz].to_vec();
        let local_vals = values[start_nnz..end_nnz].to_vec();

        SparseCsrMatrix::from_csr(end - start, matrix.ncols(), local_rp, local_ci, local_vals)
    }

    fn slice_csr_rows_owned_cols(
        matrix: &SparseCsrMatrix<S>,
        row_start: usize,
        row_end: usize,
    ) -> SparseCsrMatrix<S> {
        let row_ptr = matrix.row_ptr();
        let col_idx = matrix.col_idx();
        let values = matrix.values();
        let local_n = row_end - row_start;

        let mut local_rp = Vec::with_capacity(local_n + 1);
        let mut local_ci = Vec::new();
        let mut local_vals = Vec::new();
        local_rp.push(0);

        for global_r in row_start..row_end {
            let mut row_nnz = 0usize;
            for nz in row_ptr[global_r]..row_ptr[global_r + 1] {
                let global_c = col_idx[nz];
                if (row_start..row_end).contains(&global_c) {
                    local_ci.push(global_c - row_start);
                    local_vals.push(values[nz]);
                    row_nnz += 1;
                }
            }
            local_rp.push(local_rp.last().copied().unwrap_or(0) + row_nnz);
        }

        SparseCsrMatrix::from_csr(local_n, local_n, local_rp, local_ci, local_vals)
    }

    fn partition_rows_balanced_for_size(n_global: usize, size: usize) -> Vec<usize> {
        let base = n_global / size;
        let rem = n_global % size;
        let mut out = Vec::with_capacity(size + 1);
        out.push(0);
        for rank in 0..size {
            let take = base + usize::from(rank < rem);
            out.push(out[rank] + take);
        }
        out
    }

    fn classify_backend(size: usize, local_n: usize, global_n: usize) -> CsrBackend {
        if size <= 1 {
            CsrBackend::Serial
        } else if local_n == global_n {
            CsrBackend::Replicated
        } else {
            CsrBackend::Distributed
        }
    }

    fn try_load_rhs_s(mat_path: &Path, n: usize) -> Option<Vec<S>> {
        let mut rhs_path = mat_path.to_path_buf();
        if let Some(stem) = mat_path.file_stem().and_then(|s| s.to_str()) {
            rhs_path.set_file_name(format!("{stem}_rhs.mtx"));
            if rhs_path.exists() {
                if let Ok(mm) = read_matrix_market(&rhs_path) {
                    if let Ok(v) = mm.to_vector_scalar() {
                        if v.len() == n {
                            return Some(v);
                        }
                    }
                }
            }
        }
        None
    }

    fn assert_shifted_poisson_valid_nonzero_diagonal(matrix: &SparseCsrMatrix<S>) {
        assert_eq!(matrix.nrows(), matrix.ncols());
        assert_eq!(matrix.row_ptr().len(), matrix.nrows() + 1);
        assert!(matrix.values().len() >= matrix.nrows());

        for row in 0..matrix.nrows() {
            let start = matrix.row_ptr()[row];
            let end = matrix.row_ptr()[row + 1];
            let diag = (start..end)
                .find(|&nz| matrix.col_idx()[nz] == row)
                .map(|nz| matrix.values()[nz]);
            assert!(
                diag.is_some_and(|v| v.abs() > 0.0),
                "row {row} must have a nonzero diagonal"
            );
        }
    }

    #[test]
    fn shifted_poisson_generated_easy_and_hard_matrices_are_valid() {
        for (alpha, beta) in [(1.0, 0.25), (0.01, 0.25)] {
            let matrix = build_shifted_poisson_csr(4, 3, alpha, beta, 0.0, 0.0);
            assert_eq!(matrix.nrows(), 12);
            assert_eq!(matrix.ncols(), 12);
            assert_shifted_poisson_valid_nonzero_diagonal(&matrix);
        }
    }

    #[test]
    fn shifted_poisson_zero_convection_preserves_existing_stencil() {
        let matrix = build_shifted_poisson_csr(3, 3, 1.0, 0.25, 0.0, 0.0);
        let row = 4;
        let start = matrix.row_ptr()[row];
        let end = matrix.row_ptr()[row + 1];
        let entries = (start..end)
            .map(|nz| (matrix.col_idx()[nz], matrix.values()[nz]))
            .collect::<Vec<_>>();

        assert_eq!(entries.len(), 5);
        assert_eq!(entries[0].0, 1);
        assert_eq!(entries[0].1.real(), -1.0);
        assert_eq!(entries[1].0, 3);
        assert_eq!(entries[1].1.real(), -1.0);
        assert_eq!(entries[2].0, 4);
        assert_eq!(entries[2].1.real(), 5.0);
        assert_eq!(entries[2].1.imag(), 0.25);
        assert_eq!(entries[3].0, 5);
        assert_eq!(entries[3].1.real(), -1.0);
        assert_eq!(entries[4].0, 7);
        assert_eq!(entries[4].1.real(), -1.0);
    }

    #[test]
    fn shifted_poisson_partition_has_nonzero_local_rows_for_typical_rank_counts() {
        let matrix = build_shifted_poisson_csr(4, 4, 1.0, 0.25, 0.0, 0.0);
        for size in [1usize, 2, 3, 4, 8] {
            let part = partition_rows_balanced_for_size(matrix.nrows(), size);
            for rank in 0..size {
                let row_start = part[rank];
                let row_end = part[rank + 1];
                let local = slice_csr_rows(&matrix, row_start, row_end);
                assert!(
                    local.nrows() > 0,
                    "rank {rank}/{size} should own at least one generated Poisson row"
                );
                assert!(
                    local.values().len() > 0,
                    "rank {rank}/{size} should own generated Poisson nonzeros"
                );
            }
        }
    }

    #[test]
    fn shifted_poisson_dist_spmv_matches_global_spmv_for_ones_and_random_x() {
        #[cfg(feature = "mpi")]
        let comm = UniverseComm::Mpi(Arc::new(MpiComm::new()));
        #[cfg(not(feature = "mpi"))]
        let comm = UniverseComm::NoComm(NoComm);

        let nx = 8usize;
        let ny = 7usize;
        let alpha = 1.0;
        let beta = 0.25;
        let global_csr = build_shifted_poisson_csr(nx, ny, alpha, beta, 0.0, 0.0);
        let global_n = global_csr.nrows();
        if global_n < comm.size() {
            eprintln!(
                "skipping shifted Poisson distributed SpMV check: {global_n} rows for {} ranks",
                comm.size()
            );
            return;
        }

        let row_part = DistCsrOp::partition_rows_balanced(global_n, &comm);
        let global_row_start = row_part[comm.rank()];
        let global_row_end = row_part[comm.rank() + 1];
        let local_n = global_row_end - global_row_start;
        let local_csr = slice_csr_rows(&global_csr, global_row_start, global_row_end);
        let op = DistCsrOp::from_local_rows(
            global_n,
            global_row_start,
            &local_csr,
            &row_part,
            comm.clone(),
        )
        .expect("build distributed shifted Poisson CSR operator");

        let ones = vec![S::one(); global_n];
        let ramp = (0..global_n)
            .map(|i| S::from_parts(0.125 * (i as f64 + 1.0), -0.0625 * (i as f64 + 3.0)))
            .collect::<Vec<_>>();
        let mut state = 0x9E37_79B9_7F4A_7C15_u64;
        let random = (0..global_n)
            .map(|_| {
                state = state
                    .wrapping_mul(6364136223846793005)
                    .wrapping_add(1442695040888963407);
                let re = ((state >> 33) as f64) / ((1_u64 << 31) as f64) - 0.5;
                state = state
                    .wrapping_mul(6364136223846793005)
                    .wrapping_add(1442695040888963407);
                let im = ((state >> 33) as f64) / ((1_u64 << 31) as f64) - 0.5;
                S::from_parts(re, im)
            })
            .collect::<Vec<_>>();

        let cases = [
            ("ones", ones),
            ("complex ramp", ramp),
            ("pseudo-random", random),
        ];
        let mut scratch = BridgeScratch::default();
        for (label, x_global) in cases {
            let x_local = x_global[global_row_start..global_row_end].to_vec();
            let mut y_dist = vec![S::zero(); local_n];
            op.matvec_s(&x_local, &mut y_dist, &mut scratch);

            let mut y_global = vec![S::zero(); global_n];
            global_csr.spmv(&x_global, &mut y_global);
            let expected = &y_global[global_row_start..global_row_end];

            assert_eq!(
                y_dist.len(),
                expected.len(),
                "{label} local result length mismatch on rank {}",
                comm.rank()
            );
            for (local_row, (&actual, &want)) in y_dist.iter().zip(expected).enumerate() {
                let diff = (actual - want).abs();
                assert!(
                    diff <= 1.0e-12,
                    "{label} SpMV mismatch on rank {}, local row {local_row} (global row {}): actual={actual:?}, expected={want:?}, diff={diff:e}",
                    comm.rank(),
                    global_row_start + local_row
                );
            }
        }
    }

    #[test]
    fn shifted_poisson_problem_rhs_matches_a_times_ones() {
        #[cfg(feature = "mpi")]
        let comm = UniverseComm::Mpi(Arc::new(MpiComm::new()));
        #[cfg(not(feature = "mpi"))]
        let comm = UniverseComm::NoComm(NoComm);

        let nx = 5usize;
        let ny = 4usize;
        let alpha = 1.0;
        let beta = 0.25;
        let problem = build_shifted_poisson_case(&comm, nx, ny, alpha, beta, 0.0, 0.0)
            .expect("build generated shifted Poisson problem");
        let matrix = build_shifted_poisson_csr(nx, ny, alpha, beta, 0.0, 0.0);
        let ones = vec![S::one(); matrix.ncols()];
        let mut expected_global = vec![S::zero(); matrix.nrows()];
        matrix.spmv(&ones, &mut expected_global);
        let expected =
            &expected_global[problem.global_row_start..problem.global_row_start + problem.local_n];

        assert_eq!(problem.global_n, nx * ny);
        assert_eq!(problem.local_n, problem.rhs.len());
        assert_eq!(problem.rhs_source, RhsSource::GeneratedAOnes);
        assert!(problem.generated_case);
        assert_eq!(problem.rhs, expected);
    }

    #[test]
    fn shifted_poisson_default_pcs_include_requested_generated_set() {
        #[cfg(feature = "mpi")]
        let comm = UniverseComm::Mpi(Arc::new(MpiComm::new()));
        #[cfg(not(feature = "mpi"))]
        let comm = UniverseComm::NoComm(NoComm);

        let problem = build_shifted_poisson_case(&comm, 4, 4, 1.0, 0.25, 0.0, 0.0)
            .expect("build generated shifted Poisson problem");
        let cfg = BenchmarkConfig {
            restarts: vec![20],
            ..BenchmarkConfig::default()
        };
        let pcs = RunSpec::build_default_matrix(&cfg, &problem)
            .into_iter()
            .map(|spec| spec.pc)
            .collect::<Vec<_>>();

        assert!(pcs.contains(&PcKind::None));
        assert!(pcs.contains(&PcKind::JacobiWeak));
        assert!(pcs.contains(&PcKind::LocalIluk { k: 1 }));
        if problem.comm.size() > 1 {
            assert!(pcs.contains(&PcKind::MpiBlockJacobiIlu0Local));
        } else {
            assert!(pcs.contains(&PcKind::Ilu0Local));
        }
    }

    #[test]
    fn qc324_owned_pc_block_is_local_square_on_4_way_partition() {
        let mat_path = Path::new(env!("CARGO_MANIFEST_DIR")).join("examples/mtx/qc324.mtx");
        if !mat_path.exists() {
            eprintln!("qc324.mtx unavailable; skipping local PC block-shape check.");
            return;
        }

        let mm = read_matrix_market(&mat_path).expect("read qc324");
        let csr_sparse: SparseCsrMatrix<S> = mm.to_csr_matrix_scalar().expect("qc324 to CSR");
        let n = csr_sparse.nrows();
        let part = partition_rows_balanced_for_size(n, 4);

        for rank in 0..4 {
            let row_start = part[rank];
            let row_end = part[rank + 1];
            let local_n = row_end - row_start;
            let pc_block = slice_csr_rows_owned_cols(&csr_sparse, row_start, row_end);
            eprintln!(
                "qc324 rank {rank}: rows=[{row_start}..{row_end}), pc_block=({} x {})",
                pc_block.nrows(),
                pc_block.ncols()
            );
            assert_eq!(pc_block.nrows(), local_n, "rank {rank} local PC rows");
            assert_eq!(pc_block.ncols(), local_n, "rank {rank} local PC cols");
            assert!(
                pc_block.col_idx().iter().all(|&c| c < local_n),
                "rank {rank} has out-of-range local column index"
            );
        }
    }

    #[cfg(feature = "mpi")]
    #[test]
    fn replicated_full_ilu_apply_matches_serial_reference_under_mpi() {
        let comm = UniverseComm::Mpi(Arc::new(MpiComm::new()));
        let rank = comm.rank();
        let size = comm.size();
        let global_n = 4usize;
        let row_part = partition_rows_balanced_for_size(global_n, size);
        let row_start = row_part[rank];
        let row_end = row_part[rank + 1];
        let local_n = row_end - row_start;

        let matrix = SparseCsrMatrix::from_csr(
            global_n,
            global_n,
            vec![0, 1, 2, 3, 4],
            vec![0, 1, 2, 3],
            vec![
                S::from_parts(2.0, 0.0),
                S::from_parts(4.0, 0.0),
                S::from_parts(5.0, 0.0),
                S::from_parts(10.0, 0.0),
            ],
        );
        let x_global = vec![
            S::from_parts(2.0, 1.0),
            S::from_parts(8.0, -4.0),
            S::from_parts(15.0, 5.0),
            S::from_parts(40.0, -10.0),
        ];
        let x_local = x_global[row_start..row_end].to_vec();

        let mut cfg = IluCsrConfig::default();
        cfg.kind = IluKind::Ilu0;
        let mut serial_ilu = IluCsr::new_with_config(cfg.clone());
        serial_ilu.setup(&matrix).expect("serial ILU setup");
        let mut replicated_ilu = IluCsr::new_with_config(cfg);
        replicated_ilu.setup(&matrix).expect("replicated ILU setup");

        let mut scratch = BridgeScratch::default();
        let mut expected_global = vec![S::zero(); global_n];
        serial_ilu
            .apply_s(PcSide::Right, &x_global, &mut expected_global, &mut scratch)
            .expect("serial ILU apply");

        let mut y_local = vec![S::zero(); local_n];
        let scratch_in = std::sync::Mutex::new(vec![S::zero(); global_n]);
        let scratch_out = std::sync::Mutex::new(vec![S::zero(); global_n]);
        apply_replicated_full_ilu_owned_segment(
            &replicated_ilu,
            &comm,
            PcSide::Right,
            &x_local,
            &mut y_local,
            &mut scratch,
            global_n,
            row_start,
            local_n,
            &scratch_in,
            &scratch_out,
        )
        .expect("replicated ILU apply");

        let expected_local = &expected_global[row_start..row_end];
        assert_eq!(y_local.len(), expected_local.len());
        for (got, expected) in y_local.iter().zip(expected_local.iter()) {
            assert!(
                (*got - *expected).abs() <= 1.0e-12,
                "rank {rank} replicated apply mismatch: got={got:?}, expected={expected:?}"
            );
        }
    }

    #[test]
    fn replicated_full_ilu0_setup_helper_uses_global_dimensions_for_shifted_poisson() {
        #[cfg(feature = "mpi")]
        let comm = UniverseComm::Mpi(Arc::new(MpiComm::new()));
        #[cfg(not(feature = "mpi"))]
        let comm = UniverseComm::NoComm(NoComm);

        let problem = build_shifted_poisson_case(&comm, 5, 5, 1.0, 0.25, 0.0, 0.0)
            .expect("build generated shifted Poisson problem");
        assert_eq!(problem.global_n, 25);
        assert_eq!(
            problem.global_n % 4,
            1,
            "test matrix should not divide evenly over 4 ranks"
        );
        assert_eq!(problem.csr_for_pc.nrows(), problem.local_n);
        assert_eq!(problem.csr_for_pc.ncols(), problem.local_n);
        assert_eq!(problem.global_csr.nrows(), problem.global_n);
        assert_eq!(problem.global_csr.ncols(), problem.global_n);

        let spec = RunSpec {
            restart: 5,
            ksp: KspKind::Gmres,
            variant: FgmresVariant::Classical,
            residual_check_policy: ResidualCheckPolicy::OnConvergence,
            orthog: OrthogMethod::ClassicalGS,
            reorth: ReorthPolicy::IfNeeded,
            pc_side: PcSide::Right,
            pc: PcKind::ReplicatedFullIlu0,
        };
        let bench_cfg = BenchmarkConfig {
            run_mode: RunMode::Scalability,
            warmup_runs: 0,
            measured_runs: 0,
            ..BenchmarkConfig::default()
        };
        let pc_csr = problem.csr_for_pc.clone();
        let global_pc_csr = problem.global_csr.clone();
        let csr_pc_diag = csr_for_pc_diagnostics(
            pc_csr.as_ref(),
            problem.local_rows_nnz,
            problem.zero_global_rows_local,
        );

        let pc = setup_preconditioner_for_run_once(
            &problem,
            &spec,
            &bench_cfg,
            &pc_csr,
            &global_pc_csr,
            &csr_pc_diag,
        )
        .expect("replicated full ILU0 setup through run_once helper")
        .expect("replicated full ILU0 preconditioner handle");

        assert!(
            matches!(pc, PcHandle::ReplicatedFull { .. }),
            "expected replicated full preconditioner handle"
        );
        assert_eq!(
            KPreconditioner::dims(&pc),
            (problem.global_n, problem.global_n)
        );

        let x_local = (0..problem.local_n)
            .map(|i| S::from_parts(1.0 + (problem.global_row_start + i) as f64, -0.25))
            .collect::<Vec<_>>();
        let mut y_local = vec![S::zero(); problem.local_n];
        let mut scratch = BridgeScratch::default();
        pc.apply_s(PcSide::Right, &x_local, &mut y_local, &mut scratch)
            .expect(
                "replicated full ILU0 applies to owned local vector without dimension mismatch",
            );
    }

    #[test]
    fn replicated_full_ilu_dispatches_separately_from_local_ilu() {
        assert_eq!(
            PcKind::ReplicatedFullIlu0.dispatch_branch(),
            PcDispatchBranch::ReplicatedFullIlu
        );
        assert_eq!(
            PcKind::ReplicatedFullIluk { k: 2 }.dispatch_branch(),
            PcDispatchBranch::ReplicatedFullIlu
        );
        assert_eq!(
            PcKind::Ilu0Local.dispatch_branch(),
            PcDispatchBranch::Ilu0Local
        );
        assert_eq!(
            PcKind::MpiBlockJacobiIlu0Local.dispatch_branch(),
            PcDispatchBranch::Ilu0Local
        );
        assert_eq!(
            PcKind::LocalIluk { k: 1 }.dispatch_branch(),
            PcDispatchBranch::Ilu0Local
        );
        assert_eq!(
            PcKind::IlutLocal.dispatch_branch(),
            PcDispatchBranch::IlutLocal
        );
    }

    #[test]
    fn replicated_full_ilu_uses_global_csr_even_when_local_pc_block_is_rectangular() {
        #[cfg(feature = "mpi")]
        let comm = UniverseComm::Mpi(Arc::new(MpiComm::new()));
        #[cfg(not(feature = "mpi"))]
        let comm = UniverseComm::NoComm(NoComm);

        let op_local = SparseCsrMatrix::from_csr(
            2,
            2,
            vec![0, 1, 2],
            vec![0, 1],
            vec![S::from_parts(2.0, 0.0), S::from_parts(3.0, 0.0)],
        );
        let row_part = vec![0, 2];
        let op = DistCsrOp::from_local_rows(2, 0, &op_local, &row_part, comm.clone())
            .expect("build local dist op");

        let rectangular_pc = SparseCsrMatrix::from_csr(
            2,
            3,
            vec![0, 1, 2],
            vec![0, 1],
            vec![S::from_parts(2.0, 0.0), S::from_parts(3.0, 0.0)],
        );
        let problem = Problem {
            op: Arc::new(op),
            dist_plan_diagnostics: DistributedPlanDiagnostics {
                overlap_mode: kryst::matrix::dist_csr::HaloOverlapMode::Disabled,
                kernel_strategy: kryst::matrix::dist_csr::DistLocalKernelStrategy::RowSplitScalar,
                local_spmv_kernel: None,
                row_locality_ratio: 1.0,
                border_ratio: 0.0,
                halo_recv_volume: 0,
                halo_send_volume: 0,
                expected_communication_fraction: 0.0,
                expected_computation_fraction: 1.0,
            },
            rhs: vec![S::from_parts(2.0, 0.0), S::from_parts(3.0, 0.0)],
            rhs_source: RhsSource::GeneratedAOnes,
            solution_reference: solution_reference_diagnostics(&op_local),
            csr_for_pc: Arc::new(rectangular_pc),
            global_csr: Arc::new(op_local.clone()),
            local_rows_nnz: op_local.values().len(),
            zero_global_rows_local: count_zero_rows(&op_local),
            local_n: 2,
            global_n: 2,
            global_row_start: 0,
            comm,
            backend: CsrBackend::Serial,
            backend_descr: "unit-test".to_string(),
            generated_case: false,
        };
        let spec = RunSpec {
            restart: 5,
            ksp: KspKind::Gmres,
            variant: FgmresVariant::Classical,
            residual_check_policy: ResidualCheckPolicy::OnConvergence,
            orthog: OrthogMethod::ClassicalGS,
            reorth: ReorthPolicy::IfNeeded,
            pc_side: PcSide::Right,
            pc: PcKind::ReplicatedFullIlu0,
        };
        let bench_cfg = BenchmarkConfig {
            warmup_runs: 0,
            measured_runs: 1,
            maxits: 10,
            ..BenchmarkConfig::default()
        };

        run_once(&problem, &spec, &bench_cfg)
            .expect("replicated full ILU should set up from global CSR, not local PC CSR");
    }

    #[test]
    fn local_ilu_rejects_rectangular_owned_block() {
        #[cfg(feature = "mpi")]
        let comm = UniverseComm::Mpi(Arc::new(MpiComm::new()));
        #[cfg(not(feature = "mpi"))]
        let comm = UniverseComm::NoComm(NoComm);

        let op_local =
            SparseCsrMatrix::from_csr(2, 2, vec![0, 1, 2], vec![0, 1], vec![S::one(), S::one()]);
        let row_part = vec![0, 2];
        let op = DistCsrOp::from_local_rows(2, 0, &op_local, &row_part, comm.clone())
            .expect("build local dist op");

        let rectangular_pc =
            SparseCsrMatrix::from_csr(2, 3, vec![0, 1, 2], vec![0, 1], vec![S::one(), S::one()]);
        let problem = Problem {
            op: Arc::new(op),
            dist_plan_diagnostics: DistributedPlanDiagnostics {
                overlap_mode: kryst::matrix::dist_csr::HaloOverlapMode::Disabled,
                kernel_strategy: kryst::matrix::dist_csr::DistLocalKernelStrategy::RowSplitScalar,
                local_spmv_kernel: None,
                row_locality_ratio: 1.0,
                border_ratio: 0.0,
                halo_recv_volume: 0,
                halo_send_volume: 0,
                expected_communication_fraction: 0.0,
                expected_computation_fraction: 1.0,
            },
            rhs: vec![S::one(), S::one()],
            rhs_source: RhsSource::GeneratedAOnes,
            solution_reference: solution_reference_diagnostics(&rectangular_pc),
            csr_for_pc: Arc::new(rectangular_pc.clone()),
            global_csr: Arc::new(op_local.clone()),
            local_rows_nnz: op_local.values().len(),
            zero_global_rows_local: count_zero_rows(&op_local),
            local_n: 2,
            global_n: 2,
            global_row_start: 0,
            comm,
            backend: CsrBackend::Serial,
            backend_descr: "unit-test".to_string(),
            generated_case: false,
        };
        let spec = RunSpec {
            restart: 10,
            ksp: KspKind::Gmres,
            variant: FgmresVariant::Classical,
            residual_check_policy: ResidualCheckPolicy::OnConvergence,
            orthog: OrthogMethod::ClassicalGS,
            reorth: ReorthPolicy::IfNeeded,
            pc_side: PcSide::Right,
            pc: PcKind::Ilu0Local,
        };
        let bench_cfg = BenchmarkConfig {
            warmup_runs: 0,
            measured_runs: 1,
            ..BenchmarkConfig::default()
        };

        match run_once(&problem, &spec, &bench_cfg) {
            Err(KError::InvalidInput(msg)) => {
                assert_eq!(msg, "local ILU requires square owned block");
            }
            Err(other) => panic!("unexpected error variant: {other:?}"),
            Ok(_) => panic!("expected rectangular ILU error"),
        }
    }

    #[test]
    fn run_once_reports_requested_and_effective_policy_after_dist_policy_mutation() {
        #[cfg(feature = "mpi")]
        let comm = UniverseComm::Mpi(Arc::new(MpiComm::new()));
        #[cfg(not(feature = "mpi"))]
        let comm = UniverseComm::NoComm(NoComm);

        let op_local =
            SparseCsrMatrix::from_csr(2, 2, vec![0, 1, 2], vec![0, 1], vec![S::one(), S::one()]);
        let row_part = vec![0, 2];
        let op = DistCsrOp::from_local_rows(2, 0, &op_local, &row_part, comm.clone())
            .expect("build local dist op");

        let problem = Problem {
            op: Arc::new(op),
            dist_plan_diagnostics: DistributedPlanDiagnostics {
                overlap_mode: kryst::matrix::dist_csr::HaloOverlapMode::Disabled,
                kernel_strategy: kryst::matrix::dist_csr::DistLocalKernelStrategy::RowSplitScalar,
                local_spmv_kernel: None,
                row_locality_ratio: 1.0,
                border_ratio: 0.0,
                halo_recv_volume: 0,
                halo_send_volume: 0,
                expected_communication_fraction: 0.1,
                expected_computation_fraction: 0.5,
            },
            rhs: vec![S::one(), S::one()],
            rhs_source: RhsSource::GeneratedAOnes,
            solution_reference: solution_reference_diagnostics(&op_local),
            csr_for_pc: Arc::new(op_local.clone()),
            global_csr: Arc::new(op_local.clone()),
            local_rows_nnz: 2,
            zero_global_rows_local: count_zero_rows(&op_local),
            local_n: 2,
            global_n: 2,
            global_row_start: 0,
            comm,
            backend: CsrBackend::Serial,
            backend_descr: "unit-test".to_string(),
            generated_case: false,
        };
        let spec = RunSpec {
            restart: 5,
            ksp: KspKind::Fgmres,
            variant: FgmresVariant::Classical,
            residual_check_policy: ResidualCheckPolicy::OnConvergence,
            orthog: OrthogMethod::ClassicalGS,
            reorth: ReorthPolicy::IfNeeded,
            pc_side: PcSide::Right,
            pc: PcKind::None,
        };
        let bench_cfg = BenchmarkConfig {
            warmup_runs: 0,
            measured_runs: 1,
            ksp: KspKind::Fgmres,
            dist_policy: DistPolicyMode::Auto,
            ..BenchmarkConfig::default()
        };

        let row = run_once(&problem, &spec, &bench_cfg).expect("run once");
        assert!(row.requested_policy.contains("restart=5"));
        assert!(
            row.requested_policy
                .contains("residual-check=on-convergence")
        );
        assert!(row.effective_policy.contains("restart=16"));
        assert!(
            row.effective_policy
                .contains("residual-check=every-iteration")
        );

        let rendered = render_result_row(&row, RunMode::Correctness, &problem);
        assert!(rendered.contains(&row.requested_policy));
        assert!(rendered.contains(&row.effective_policy));
        assert_eq!(row.pc_apply, "none");
        assert_eq!(row.dist_ok, Some(true));
        assert_eq!(row.global_ok, Some(true));
        assert_eq!(row.x_ok, Some(true));
        assert!(rendered.contains("OK       OK       OK"));
        assert!(rendered.contains("csr-cx"));
        assert!(rendered.contains("ser"));
        assert!(rendered.contains("none"));
    }

    #[test]
    fn generated_aones_empty_row_suppresses_x_error_rel() {
        #[cfg(feature = "mpi")]
        let comm = UniverseComm::Mpi(Arc::new(MpiComm::new()));
        #[cfg(not(feature = "mpi"))]
        let comm = UniverseComm::NoComm(NoComm);

        let op_local = SparseCsrMatrix::from_csr(2, 2, vec![0, 1, 1], vec![0], vec![S::one()]);
        let row_part = vec![0, 2];
        let op = DistCsrOp::from_local_rows(2, 0, &op_local, &row_part, comm.clone())
            .expect("build local dist op");
        let diagnostics = solution_reference_diagnostics(&op_local);
        assert_eq!(diagnostics.zero_row_count, 1);
        assert_eq!(diagnostics.missing_diagonal_count, 1);
        assert!(!diagnostics.valid_for_x_error());

        let problem = Problem {
            op: Arc::new(op),
            dist_plan_diagnostics: DistributedPlanDiagnostics {
                overlap_mode: kryst::matrix::dist_csr::HaloOverlapMode::Disabled,
                kernel_strategy: kryst::matrix::dist_csr::DistLocalKernelStrategy::RowSplitScalar,
                local_spmv_kernel: None,
                row_locality_ratio: 1.0,
                border_ratio: 0.0,
                halo_recv_volume: 0,
                halo_send_volume: 0,
                expected_communication_fraction: 0.0,
                expected_computation_fraction: 1.0,
            },
            rhs: vec![S::one(), S::zero()],
            rhs_source: RhsSource::GeneratedAOnes,
            solution_reference: diagnostics,
            csr_for_pc: Arc::new(op_local.clone()),
            global_csr: Arc::new(op_local.clone()),
            local_rows_nnz: 1,
            zero_global_rows_local: count_zero_rows(&op_local),
            local_n: 2,
            global_n: 2,
            global_row_start: 0,
            comm,
            backend: CsrBackend::Serial,
            backend_descr: "unit-test".to_string(),
            generated_case: false,
        };
        let spec = RunSpec {
            restart: 5,
            ksp: KspKind::Gmres,
            variant: FgmresVariant::Classical,
            residual_check_policy: ResidualCheckPolicy::OnConvergence,
            orthog: OrthogMethod::ClassicalGS,
            reorth: ReorthPolicy::IfNeeded,
            pc_side: PcSide::Right,
            pc: PcKind::None,
        };
        let bench_cfg = BenchmarkConfig {
            run_mode: RunMode::Correctness,
            warmup_runs: 0,
            measured_runs: 1,
            ..BenchmarkConfig::default()
        };

        let row = run_once(&problem, &spec, &bench_cfg).expect("run once");
        assert_eq!(problem.rhs_source, RhsSource::GeneratedAOnes);
        assert_eq!(row.x_error_rel, None);
        assert_eq!(row.global_x_error_rel, None);
        assert_eq!(row.x_ok, None);
        let rendered = render_result_row(&row, RunMode::Correctness, &problem);
        assert!(rendered.contains("OK       OK      N/A"));
    }

    #[test]
    fn correctness_default_dist_policy_off_preserves_requested_restart() {
        #[cfg(feature = "mpi")]
        let comm = UniverseComm::Mpi(Arc::new(MpiComm::new()));
        #[cfg(not(feature = "mpi"))]
        let comm = UniverseComm::NoComm(NoComm);

        let op_local =
            SparseCsrMatrix::from_csr(2, 2, vec![0, 1, 2], vec![0, 1], vec![S::one(), S::one()]);
        let row_part = vec![0, 2];
        let op = DistCsrOp::from_local_rows(2, 0, &op_local, &row_part, comm.clone())
            .expect("build local dist op");

        let problem = Problem {
            op: Arc::new(op),
            dist_plan_diagnostics: DistributedPlanDiagnostics {
                overlap_mode: kryst::matrix::dist_csr::HaloOverlapMode::Disabled,
                kernel_strategy: kryst::matrix::dist_csr::DistLocalKernelStrategy::RowSplitScalar,
                local_spmv_kernel: None,
                row_locality_ratio: 1.0,
                border_ratio: 0.0,
                halo_recv_volume: 0,
                halo_send_volume: 0,
                expected_communication_fraction: 0.1,
                expected_computation_fraction: 0.5,
            },
            rhs: vec![S::one(), S::one()],
            rhs_source: RhsSource::GeneratedAOnes,
            solution_reference: solution_reference_diagnostics(&op_local),
            csr_for_pc: Arc::new(op_local.clone()),
            global_csr: Arc::new(op_local.clone()),
            local_rows_nnz: 2,
            zero_global_rows_local: count_zero_rows(&op_local),
            local_n: 2,
            global_n: 2,
            global_row_start: 0,
            comm,
            backend: CsrBackend::Serial,
            backend_descr: "unit-test".to_string(),
            generated_case: false,
        };
        let spec = RunSpec {
            restart: 5,
            ksp: KspKind::Gmres,
            variant: FgmresVariant::Classical,
            residual_check_policy: ResidualCheckPolicy::OnConvergence,
            orthog: OrthogMethod::ClassicalGS,
            reorth: ReorthPolicy::IfNeeded,
            pc_side: PcSide::Right,
            pc: PcKind::None,
        };
        let bench_cfg = BenchmarkConfig {
            run_mode: RunMode::Correctness,
            warmup_runs: 0,
            measured_runs: 1,
            ..BenchmarkConfig::default()
        };

        assert_eq!(bench_cfg.dist_policy, DistPolicyMode::Off);
        let row = run_once(&problem, &spec, &bench_cfg).expect("run once");
        assert!(row.requested_policy.contains("restart=5"));
        assert!(row.effective_policy.contains("restart=5"));
        assert!(
            row.effective_policy
                .contains("residual-check=on-convergence")
        );
    }

    #[test]
    fn metadata_labels_snapshot_serial() {
        assert_eq!(
            pc_domain_label_for_size(PcKind::ReplicatedFullIlu0, 1),
            "full"
        );
        assert_eq!(
            pc_domain_label_for_size(PcKind::ReplicatedFullIluk { k: 1 }, 1),
            "full"
        );
        assert_eq!(pc_domain_label_for_size(PcKind::Ilu0Local, 1), "full");
        assert_eq!(pc_domain_label_for_size(PcKind::IlutLocal, 1), "full");
        assert_eq!(
            pc_domain_label_for_size(PcKind::LocalIluk { k: 1 }, 1),
            "full"
        );
        assert_eq!(
            pc_domain_label_for_size(PcKind::MpiBlockJacobiIlu0Local, 1),
            "own0"
        );
        assert_eq!(pc_domain_label_for_size(PcKind::None, 1), "n/a");
        assert_eq!(pc_domain_label_for_size(PcKind::JacobiWeak, 1), "n/a");
        assert_eq!(
            pc_apply_label_for_size(PcKind::ReplicatedFullIlu0, 1),
            "replicated/full"
        );
    }

    #[test]
    fn metadata_labels_snapshot_mpi() {
        assert_eq!(
            pc_domain_label_for_size(PcKind::ReplicatedFullIlu0, 2),
            "full"
        );
        assert_eq!(
            pc_apply_label_for_size(PcKind::ReplicatedFullIlu0, 2),
            "replicated/full/allgather, not scalable"
        );
        assert_eq!(
            pc_domain_label_for_size(PcKind::ReplicatedFullIluk { k: 1 }, 2),
            "full"
        );
        assert_eq!(
            pc_apply_label_for_size(PcKind::ReplicatedFullIluk { k: 1 }, 2),
            "replicated/full/allgather, not scalable"
        );
        assert_eq!(pc_domain_label_for_size(PcKind::Ilu0Local, 2), "own0");
        assert_eq!(pc_domain_label_for_size(PcKind::IlutLocal, 2), "own0");
        assert_eq!(
            pc_domain_label_for_size(PcKind::LocalIluk { k: 1 }, 2),
            "own0"
        );
        assert_eq!(
            pc_domain_label_for_size(PcKind::MpiBlockJacobiIlu0Local, 2),
            "own0"
        );
        assert_eq!(
            pc_apply_label_for_size(PcKind::Ilu0Local, 2),
            "owned-block/overlap0"
        );
        assert_eq!(
            pc_domain_label_for_size(PcKind::AsmIlu0Overlap { overlap: 1 }, 2),
            "own+ghost1"
        );
        assert_eq!(
            pc_domain_label_for_size(PcKind::AsmIlu0Overlap { overlap: 2 }, 2),
            "own+ghost2"
        );
        assert_eq!(
            pc_apply_label_for_size(PcKind::AsmIlu0Overlap { overlap: 1 }, 2),
            "asm/overlap-gather/owned-output"
        );
        assert_eq!(
            pc_apply_label_for_size(PcKind::RasIlu0Overlap { overlap: 1 }, 2),
            "ras/overlap-gather/restricted-owned"
        );
        assert_eq!(pc_domain_label_for_size(PcKind::None, 2), "n/a");
        assert_eq!(pc_apply_label_for_size(PcKind::None, 2), "none");
        assert_eq!(
            pc_apply_label_for_size(PcKind::JacobiWeak, 2),
            "jacobi/local"
        );
    }

    #[test]
    fn row_scaled_mode_scales_pc_setup_matrix_rows() {
        let pc = SparseCsrMatrix::from_csr(
            2,
            2,
            vec![0, 2, 4],
            vec![0, 1, 0, 1],
            vec![
                S::new(2.0, 0.0),
                S::new(0.0, 0.0),
                S::new(0.0, 0.0),
                S::new(4.0, 0.0),
            ],
        );
        let row_scale = compute_row_scaling(&pc, 1e-15);
        let pc_scaled = scale_csr_rows(&pc, &row_scale);

        assert_eq!(row_scale, vec![0.5, 0.25]);
        assert_eq!(pc_scaled.values()[0], S::one());
        assert_eq!(pc_scaled.values()[3], S::one());
        assert_ne!(pc_scaled.values(), pc.values());
    }

    #[test]
    fn mpi_mode_run_matrix_has_no_duplicate_semantic_pc_experiments() {
        let cfg = BenchmarkConfig::default();
        let mpi_mode = true;
        let pcs = if cfg.pcs.is_empty() {
            match cfg.run_mode {
                RunMode::Correctness => {
                    if mpi_mode {
                        vec![
                            PcKind::ReplicatedFullIlu0,
                            PcKind::None,
                            PcKind::JacobiWeak,
                            PcKind::MpiBlockJacobiIlu0Local,
                        ]
                    } else {
                        vec![PcKind::Ilu0Local, PcKind::JacobiWeak, PcKind::None]
                    }
                }
                RunMode::Scalability => {
                    if mpi_mode {
                        vec![PcKind::MpiBlockJacobiIlu0Local]
                    } else {
                        vec![PcKind::Ilu0Local]
                    }
                }
                RunMode::Robustness => {
                    if mpi_mode {
                        vec![PcKind::MpiBlockJacobiIlu0Local, PcKind::JacobiWeak]
                    } else {
                        vec![PcKind::Ilu0Local, PcKind::JacobiWeak]
                    }
                }
            }
        } else {
            cfg.pcs.clone()
        };

        assert_eq!(
            pcs,
            vec![
                PcKind::ReplicatedFullIlu0,
                PcKind::None,
                PcKind::JacobiWeak,
                PcKind::MpiBlockJacobiIlu0Local,
            ]
        );

        let mut seen = std::collections::BTreeSet::new();
        for &restart in &cfg.restarts {
            for &variant in &cfg.variants {
                for &orthog in &cfg.orthogs {
                    for &reorth in &cfg.reorths {
                        for &pc in &pcs {
                            let key = format!(
                                "restart={restart}|variant={}|orthog={}|reorth={}|pc={}",
                                variant_label(variant),
                                orthog_label(orthog),
                                reorth_label(reorth),
                                pc.semantic_experiment_key(mpi_mode)
                            );
                            assert!(
                                seen.insert(key.clone()),
                                "duplicate semantic benchmark row: {key}"
                            );
                        }
                    }
                }
            }
        }
    }

    #[test]
    fn parse_pc_accepts_documented_sweep_names() {
        assert_eq!(
            parse_pc_csv(
                "--pcs",
                "replicated-full-ilu0,none,jacobi,block-jacobi-ilu0"
            )
            .expect("parse documented correctness preconditioners"),
            vec![
                PcKind::ReplicatedFullIlu0,
                PcKind::None,
                PcKind::JacobiWeak,
                PcKind::MpiBlockJacobiIlu0Local,
            ]
        );
        assert_eq!(
            parse_pc_csv("--pcs", "none,jacobi,block-jacobi-ilu0")
                .expect("parse documented scalability preconditioners"),
            vec![
                PcKind::None,
                PcKind::JacobiWeak,
                PcKind::MpiBlockJacobiIlu0Local,
            ]
        );
    }

    #[test]
    fn parse_pc_accepts_replicated_full_ilu0_alias() {
        assert_eq!(
            parse_pc("replicated-full-ilu0").expect("parse replicated-full-ilu0"),
            PcKind::ReplicatedFullIlu0
        );
        assert_eq!(
            parse_pc("replicated-ilu0").expect("parse replicated-ilu0"),
            PcKind::ReplicatedFullIlu0
        );
    }

    #[test]
    fn parse_pc_accepts_block_jacobi_iluk_aliases() {
        assert_eq!(
            parse_pc("block-jacobi-iluk1").expect("parse block-jacobi-iluk1"),
            PcKind::LocalIluk { k: 1 }
        );
        assert_eq!(
            parse_pc("block-jacobi-iluk1-overlap0").expect("parse block-jacobi-iluk1-overlap0"),
            PcKind::LocalIluk { k: 1 }
        );
        assert_eq!(
            parse_pc("block-jacobi-iluk:2").expect("parse block-jacobi-iluk:2"),
            PcKind::LocalIluk { k: 2 }
        );
        assert_eq!(
            parse_pc("block-jacobi-iluk3-overlap0").expect("parse block-jacobi-iluk3-overlap0"),
            PcKind::LocalIluk { k: 3 }
        );
    }

    #[test]
    fn parse_pc_accepts_overlap_aliases() {
        assert_eq!(
            parse_pc("asm-ilu0-overlap1").expect("parse asm-ilu0-overlap1"),
            PcKind::AsmIlu0Overlap { overlap: 1 }
        );
        assert_eq!(
            parse_pc("asm-ilu0-overlap2").expect("parse asm-ilu0-overlap2"),
            PcKind::AsmIlu0Overlap { overlap: 2 }
        );
        assert_eq!(
            parse_pc("ras-ilu0-overlap1").expect("parse ras-ilu0-overlap1"),
            PcKind::RasIlu0Overlap { overlap: 1 }
        );
        assert_eq!(
            parse_pc("ras-iluk1-overlap1").expect("parse ras-iluk1-overlap1"),
            PcKind::RasIlukOverlap { k: 1, overlap: 1 }
        );
    }

    #[test]
    fn local_iluk_label_renders_block_jacobi_alias() {
        assert_eq!(
            PcKind::LocalIluk { k: 3 }.label(),
            "block-jacobi-iluk3-overlap0"
        );
    }

    #[test]
    fn preconditioner_dispatch_variants_are_distinct_or_explicit_aliases() {
        let variants = [
            PcKind::None,
            PcKind::JacobiWeak,
            PcKind::Ilu0Local,
            PcKind::IlutLocal,
            PcKind::MpiBlockJacobiIlu0Local,
        ];

        for (idx, left) in variants.iter().copied().enumerate() {
            for right in variants.iter().copied().skip(idx + 1) {
                if left.dispatch_branch() == right.dispatch_branch() {
                    let explicit_alias = left.explicit_alias_of() == Some(right)
                        || right.explicit_alias_of() == Some(left);
                    assert!(
                        explicit_alias,
                        "dispatch collision without alias marker: {left:?} vs {right:?}"
                    );
                }
            }
        }
    }

    #[test]
    fn method_label_distinguishes_gmres_and_fgmres() {
        let base = RunSpec {
            restart: 50,
            ksp: KspKind::Gmres,
            variant: FgmresVariant::Classical,
            residual_check_policy: ResidualCheckPolicy::OnConvergence,
            orthog: OrthogMethod::ClassicalGS,
            reorth: ReorthPolicy::IfNeeded,
            pc_side: PcSide::Right,
            pc: PcKind::JacobiWeak,
        };
        assert!(base.method_label().starts_with("GMRES+"));

        let flexible = RunSpec {
            ksp: KspKind::Fgmres,
            ..base
        };
        assert!(flexible.method_label().starts_with("FGMRES+"));
    }

    #[test]
    fn fixed_pc_default_run_specs_select_gmres() {
        #[cfg(feature = "mpi")]
        let comm = UniverseComm::Mpi(Arc::new(MpiComm::new()));
        #[cfg(not(feature = "mpi"))]
        let comm = UniverseComm::NoComm(NoComm);

        let problem = build_shifted_poisson_case(&comm, 4, 4, 1.0, 0.25, 0.0, 0.0)
            .expect("build generated shifted Poisson problem");
        let cfg = BenchmarkConfig {
            restarts: vec![8],
            pcs: vec![PcKind::None, PcKind::JacobiWeak, PcKind::Ilu0Local],
            ..BenchmarkConfig::default()
        };
        let runs = RunSpec::build_default_matrix(&cfg, &problem);
        assert_eq!(runs.len(), 3);
        assert!(runs.iter().all(|spec| spec.ksp == KspKind::Gmres));
    }

    #[test]
    fn ilut_complex_label_explicitly_marks_real_projection_fallback_path() {
        let label = PcKind::IlutLocal.label();
        assert!(label.contains(ILUT_REAL_PROJECTION_FALLBACK_LABEL));
        assert!(label.contains("degraded/provisional"));
        assert!(label.contains("real-projection fallback"));
    }

    #[test]
    fn ilut_method_row_and_export_labels_are_self_contained() {
        let spec = RunSpec {
            restart: 50,
            ksp: KspKind::Gmres,
            variant: FgmresVariant::Classical,
            residual_check_policy: ResidualCheckPolicy::OnConvergence,
            orthog: OrthogMethod::ClassicalGS,
            reorth: ReorthPolicy::IfNeeded,
            pc_side: PcSide::Right,
            pc: PcKind::IlutLocal,
        };

        assert!(
            spec.method_label()
                .contains(ILUT_REAL_PROJECTION_FALLBACK_LABEL)
        );
        assert!(
            PcKind::IlutLocal
                .semantic_experiment_key(false)
                .contains("real-projection-fallback")
        );
    }

    #[test]
    fn correctness_mode_disables_stagnation_fallback_and_keeps_restart_target() {
        let spec = RunSpec {
            restart: 50,
            ksp: KspKind::Gmres,
            variant: FgmresVariant::Classical,
            residual_check_policy: ResidualCheckPolicy::OnConvergence,
            orthog: OrthogMethod::ClassicalGS,
            reorth: ReorthPolicy::IfNeeded,
            pc_side: PcSide::Right,
            pc: PcKind::None,
        };
        let cfg = BenchmarkConfig {
            run_mode: RunMode::Correctness,
            allow_stagnation_fallback: true,
            min_inner_before_fallback: 12,
            ..BenchmarkConfig::default()
        };
        let solver = configured_fgmres_solver(&spec, &cfg);
        assert_eq!(solver.restart, 50);
        assert_eq!(solver.stagnation_policy, FgmresStagnationPolicy::Disabled);
        assert_eq!(solver.min_inner_before_fallback, 12);
    }

    #[test]
    fn benchmark_config_defaults_dist_policy_from_run_mode() {
        let default_cfg = BenchmarkConfig::default();
        assert_eq!(default_cfg.dist_policy, DistPolicyMode::Off);

        let scalability_cfg =
            BenchmarkConfig::from_args(vec!["--mode".to_string(), "scalability".to_string()])
                .expect("parse scalability mode");
        assert_eq!(scalability_cfg.dist_policy, DistPolicyMode::Auto);

        let robustness_cfg =
            BenchmarkConfig::from_args(vec!["--mode".to_string(), "robustness".to_string()])
                .expect("parse robustness mode");
        assert_eq!(robustness_cfg.dist_policy, DistPolicyMode::Off);
    }

    #[test]
    fn benchmark_config_parses_ksp_selection() {
        let gmres_cfg = BenchmarkConfig::from_args(vec!["--ksp".to_string(), "gmres".to_string()])
            .expect("parse --ksp gmres");
        assert_eq!(gmres_cfg.ksp, KspKind::Gmres);

        let fgmres_cfg =
            BenchmarkConfig::from_args(vec!["--ksp".to_string(), "fgmres".to_string()])
                .expect("parse --ksp fgmres");
        assert_eq!(fgmres_cfg.ksp, KspKind::Fgmres);
    }

    #[test]
    fn benchmark_config_rejects_invalid_ksp_selection() {
        let err = BenchmarkConfig::from_args(vec!["--ksp".to_string(), "bicgstab".to_string()])
            .expect_err("invalid --ksp should be rejected");
        match err {
            KError::InvalidInput(msg) => assert!(msg.contains("--ksp") || msg.contains("KSP")),
            other => panic!("unexpected error variant: {other:?}"),
        }
    }

    #[test]
    fn benchmark_config_parses_explicit_dist_policy_override() {
        let cfg = BenchmarkConfig::from_args(vec![
            "--mode".to_string(),
            "scalability".to_string(),
            "--dist-policy".to_string(),
            "off".to_string(),
        ])
        .expect("parse explicit dist policy");
        assert_eq!(cfg.dist_policy, DistPolicyMode::Off);
    }

    #[test]
    fn benchmark_config_parses_poisson_shift_and_convection_flags() {
        let cfg = BenchmarkConfig::from_args(vec![
            "--poisson-shift-real".to_string(),
            "0.01".to_string(),
            "--poisson-shift-imag".to_string(),
            "-0.5".to_string(),
            "--poisson-convection".to_string(),
            "1.25,-2.5".to_string(),
        ])
        .expect("parse shifted Poisson overrides");

        assert_eq!(cfg.poisson_shift_real, 0.01);
        assert_eq!(cfg.poisson_shift_imag, -0.5);
        assert_eq!(cfg.poisson_convection_x, 1.25);
        assert_eq!(cfg.poisson_convection_y, -2.5);
    }

    #[test]
    fn benchmark_config_rejects_invalid_poisson_convection_flag() {
        let err = BenchmarkConfig::from_args(vec![
            "--poisson-convection".to_string(),
            "1.0;2.0".to_string(),
        ])
        .expect_err("invalid convection pair should be rejected");
        match err {
            KError::InvalidInput(msg) => assert!(msg.contains("--poisson-convection")),
            other => panic!("unexpected error variant: {other:?}"),
        }
    }

    #[test]
    fn benchmark_config_parses_block_jacobi_iluk_overlap_pc() {
        let cfg = BenchmarkConfig::from_args(vec![
            "--pcs".to_string(),
            "block-jacobi-iluk1-overlap0".to_string(),
        ])
        .expect("parse block Jacobi ILU(k) pc alias");

        assert_eq!(cfg.pcs, vec![PcKind::LocalIluk { k: 1 }]);
    }

    #[test]
    fn benchmark_config_uses_conservative_default_fgmres_haptol() {
        let cfg = BenchmarkConfig::default();
        assert_eq!(cfg.fgmres_haptol, 1e-30);
    }

    #[test]
    fn benchmark_config_parses_fgmres_haptol_override() {
        let cfg =
            BenchmarkConfig::from_args(vec!["--fgmres-haptol".to_string(), "1e-22".to_string()])
                .expect("parse args");
        assert_eq!(cfg.fgmres_haptol, 1e-22);
    }

    #[test]
    fn benchmark_config_rejects_non_positive_or_non_finite_fgmres_haptol() {
        for bad in ["0", "-1", "NaN", "inf"] {
            let err =
                BenchmarkConfig::from_args(vec!["--fgmres-haptol".to_string(), bad.to_string()])
                    .expect_err("invalid --fgmres-haptol should be rejected");
            match err {
                KError::InvalidInput(msg) => assert!(msg.contains("--fgmres-haptol")),
                other => panic!("unexpected error variant: {other:?}"),
            }
        }
    }

    #[test]
    fn benchmark_config_accepts_include_overlap_pcs_flag() {
        let cfg = BenchmarkConfig::from_args(vec!["--include-overlap-pcs".to_string()])
            .expect("parse args");
        assert!(cfg.include_overlap_pcs);
    }

    #[test]
    fn include_overlap_pcs_is_limited_to_mpi_correctness_mode() {
        let correctness_cfg = BenchmarkConfig {
            include_overlap_pcs: true,
            run_mode: RunMode::Correctness,
            ..BenchmarkConfig::default()
        };
        assert!(should_include_overlap_pcs(&correctness_cfg, 4));
        assert!(!should_include_overlap_pcs(&correctness_cfg, 1));

        let scalability_cfg = BenchmarkConfig {
            include_overlap_pcs: true,
            run_mode: RunMode::Scalability,
            ..BenchmarkConfig::default()
        };
        assert!(!should_include_overlap_pcs(&scalability_cfg, 4));
    }

    #[test]
    fn include_overlap_pcs_appends_verdict_gated_mpi_rows_after_block_jacobi_baseline() {
        let mut pcs = vec![PcKind::ReplicatedFullIlu0, PcKind::None, PcKind::JacobiWeak];
        append_mpi_correctness_overlap_pcs(&mut pcs);
        assert_eq!(
            pcs,
            vec![
                PcKind::ReplicatedFullIlu0,
                PcKind::None,
                PcKind::JacobiWeak,
                PcKind::MpiBlockJacobiIlu0Local,
                PcKind::AsmIlu0Overlap { overlap: 1 },
                PcKind::RasIlu0Overlap { overlap: 1 },
                PcKind::RasIlukOverlap { k: 1, overlap: 1 },
            ]
        );
    }

    #[test]
    fn overlap_pcs_dispatch_through_overlap_ilu_branch() {
        for pc in MPI_CORRECTNESS_OVERLAP_PCS {
            assert_eq!(pc.dispatch_branch(), PcDispatchBranch::OverlapIlu);
        }
    }

    #[test]
    fn benchmark_config_accepts_mark_replicated_check_flag() {
        let cfg = BenchmarkConfig::from_args(vec!["--mark-replicated-check".to_string()])
            .expect("parse args");
        assert!(cfg.mark_replicated_check);
    }

    #[test]
    fn benchmark_config_rejects_deprecated_correctness_replicated_check_flag() {
        let err = BenchmarkConfig::from_args(vec!["--correctness-replicated-check".to_string()])
            .expect_err("deprecated flag should be rejected");
        match err {
            KError::InvalidInput(msg) => {
                assert!(msg.contains("renamed to --mark-replicated-check"))
            }
            other => panic!("unexpected error variant: {other:?}"),
        }
    }

    #[test]
    fn diagnostics_format_includes_global_and_rank_spread_fields() {
        let diag = CsrForPcDiagnosticsGlobal {
            nnz_local_block: rank_spread_from_locals(&[12.0, 30.0], 12.0),
            nnz_local_rows: rank_spread_from_locals(&[10.0, 20.0], 10.0),
            nnz_ratio: rank_spread_from_locals(&[1.2, 1.5], 1.2),
            zero_local_rows_rank_count: 1.0,
            zero_global_rows: rank_spread_from_locals(&[0.0, 2.0], 0.0),
            structural_diag_missing: rank_spread_from_locals(&[1.0, 3.0], 1.0),
            diag_min_abs: rank_spread_from_locals(&[1e-8, 1e-10], 1e-8),
            diag_max_abs: rank_spread_from_locals(&[9.0, 11.0], 9.0),
            diag_tiny_or_missing_count: rank_spread_from_locals(&[1.0, 4.0], 1.0),
            pivot_perturbation_count: Some(rank_spread_from_locals(&[0.0, 3.0], 0.0)),
        };
        let text = format_csr_for_pc_diagnostics("ilu0-local", &diag);
        assert!(text.contains("global_sum="));
        assert!(text.contains("rank_min="));
        assert!(text.contains("rank_max="));
        assert!(text.contains("rank0_local="));
        assert!(text.contains("nnz_ratio_global=1.400e0"));
        assert!(text.contains("rank_local_min=1.200e0"));
        assert!(text.contains("rank_local_max=1.500e0"));
        assert!(!text.contains("nnz_ratio:global_sum="));
        assert!(text.contains("zero_local_rows_rank_count=1"));
        assert!(text.contains("zero_global_rows:global_sum=2"));
        assert!(text.contains("structural_diag_missing:global_sum=4"));
    }

    #[test]
    fn diagnostics_capture_nonzero_off_rank_contribution() {
        let spread = rank_spread_from_locals(&[5.0, 0.0, 7.0], 5.0);
        assert_eq!(spread.rank0_local, 5.0);
        assert_eq!(spread.global_sum, 12.0);
        assert!(spread.global_sum > spread.rank0_local);
        assert_eq!(spread.rank_max, 7.0);
    }
}

#[cfg(feature = "complex")]
fn main() -> Result<(), KError> {
    complex_demo::run()
}