kryst 4.0.3

//! Convergence tracking & tolerance checks for iterative solvers.

#[allow(unused_imports)]
use crate::algebra::prelude::*;
use crate::error::KError;

/// Stable reason categories for automation/diagnostics.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum ReasonCategory {
    Breakdown,
    Nan,
    Inf,
    PcSetup,
    PcApply,
}

/// Canonical failure mapping keys used by KSP/PC code paths.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum FailureReasonKind {
    Breakdown,
    BreakdownBiCG,
    IndefiniteMatrix,
    IndefinitePc,
    Nan,
    Inf,
    PcSetup,
    PcApply,
}

/// Stage for mapping transport errors into stable KSP reasons.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum FailureStage {
    Setup,
    Solve,
}

/// Convergence criteria for iterative solvers.
///
/// This struct defines four types of stopping criteria:
/// - **Relative tolerance**: `‖r‖/‖b‖ ≤ rtol`
/// - **Absolute tolerance**: `‖r‖ ≤ atol`
/// - **Divergence threshold**: `‖r‖ ≥ dtol * ‖b‖`
/// - **Maximum iterations**: `iterations ≥ max_iters`
pub struct Convergence {
    /// Relative tolerance: ‖r‖/‖b‖ ≤ rtol ⇒ converge
    pub rtol: R,
    /// Absolute tolerance: ‖r‖ ≤ atol ⇒ converge
    pub atol: R,
    /// Divergence threshold: ‖r‖ ≥ dtol * ‖b‖ ⇒ diverge
    pub dtol: R,
    /// Maximum iterations
    pub max_iters: usize,
}

/// Reason for convergence or divergence.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum ConvergedReason {
    /// Converged due to relative tolerance: ‖r‖/‖b‖ ≤ rtol
    ConvergedRtol,
    /// Converged due to absolute tolerance: ‖r‖ ≤ atol
    ConvergedAtol,
    /// Converged because the step reached a trust-region bound
    ConvergedTrustRegion,
    /// Converged due to a happy breakdown (e.g., `pᵀAp` ≈ 0)
    ConvergedHappyBreakdown,
    /// Diverged due to NaN residuals
    DivergedNan,
    /// Diverged due to Inf residuals
    DivergedInf,
    /// Diverged due to divergence tolerance: ‖r‖ ≥ dtol * ‖b‖
    DivergedDtol,
    /// Diverged due to maximum iterations reached
    DivergedMaxIts,
    /// Diverged due to a breakdown (generic)
    DivergedBreakdown,
    /// Diverged because the Arnoldi basis lost rank during orthogonalization.
    DivergedArnoldiRankLoss,
    /// Diverged due to a BiCG-specific breakdown
    DivergedBreakdownBiCG,
    /// Diverged because the reduced Hessenberg system became singular/near-singular.
    DivergedReducedSystemSingular,
    /// Diverged because the matrix is indefinite
    DivergedIndefiniteMatrix,
    /// Diverged because the preconditioner is indefinite
    DivergedIndefinitePC,
    /// Diverged because preconditioner setup failed
    DivergedPcSetupFailed,
    /// Diverged because the preconditioner failed
    DivergedPcFailed,
    /// Diverged due to a monitor-requested stop
    StoppedByMonitor,
    /// Continue iterating (none of the stopping criteria met)
    Continued,
}

/// Contract acceptance status derived from solver reason and true residual.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum AcceptanceStatus {
    Ok,
    OkWithWarning,
    ContractMismatch,
    Breakdown,
    Stagnated,
    Failed,
}

impl AcceptanceStatus {
    pub fn as_str(self) -> &'static str {
        match self {
            AcceptanceStatus::Ok => "ok",
            AcceptanceStatus::OkWithWarning => "ok_with_warning",
            AcceptanceStatus::ContractMismatch => "contract_mismatch",
            AcceptanceStatus::Breakdown => "breakdown",
            AcceptanceStatus::Stagnated => "stagnated",
            AcceptanceStatus::Failed => "failed",
        }
    }

    pub fn is_accepted(self) -> bool {
        matches!(self, AcceptanceStatus::Ok | AcceptanceStatus::OkWithWarning)
    }
}

/// Classify acceptance from a canonical solver reason plus true residual contract.
///
/// The true residual check is evaluated first. If `true_residual < tol`, the solve is
/// accepted even if the reported reason is non-converged, but marked as `OkWithWarning`.
pub fn classify_acceptance_status(
    reason: ConvergedReason,
    true_residual: f64,
    tol: f64,
) -> AcceptanceStatus {
    let meets_true_residual = true_residual.is_finite() && true_residual < tol;
    if meets_true_residual {
        if reason.is_converged() {
            return AcceptanceStatus::Ok;
        }
        return AcceptanceStatus::OkWithWarning;
    }

    match reason {
        ConvergedReason::ConvergedRtol
        | ConvergedReason::ConvergedAtol
        | ConvergedReason::ConvergedTrustRegion
        | ConvergedReason::ConvergedHappyBreakdown
        | ConvergedReason::DivergedIndefiniteMatrix
        | ConvergedReason::DivergedIndefinitePC => AcceptanceStatus::ContractMismatch,
        ConvergedReason::DivergedBreakdown
        | ConvergedReason::DivergedArnoldiRankLoss
        | ConvergedReason::DivergedBreakdownBiCG
        | ConvergedReason::DivergedReducedSystemSingular
        | ConvergedReason::DivergedNan
        | ConvergedReason::DivergedInf => AcceptanceStatus::Breakdown,
        ConvergedReason::DivergedDtol
        | ConvergedReason::DivergedMaxIts
        | ConvergedReason::StoppedByMonitor => AcceptanceStatus::Stagnated,
        ConvergedReason::DivergedPcSetupFailed
        | ConvergedReason::DivergedPcFailed
        | ConvergedReason::Continued => AcceptanceStatus::Failed,
    }
}

impl ConvergedReason {
    /// Canonicalize a method-local failure kind into a stable reason code.
    pub fn from_failure_kind(kind: FailureReasonKind) -> Self {
        match kind {
            FailureReasonKind::Breakdown => ConvergedReason::DivergedBreakdown,
            FailureReasonKind::BreakdownBiCG => ConvergedReason::DivergedBreakdownBiCG,
            FailureReasonKind::IndefiniteMatrix => ConvergedReason::DivergedIndefiniteMatrix,
            FailureReasonKind::IndefinitePc => ConvergedReason::DivergedIndefinitePC,
            FailureReasonKind::Nan => ConvergedReason::DivergedNan,
            FailureReasonKind::Inf => ConvergedReason::DivergedInf,
            FailureReasonKind::PcSetup => ConvergedReason::DivergedPcSetupFailed,
            FailureReasonKind::PcApply => ConvergedReason::DivergedPcFailed,
        }
    }

    /// Group reason codes into stable automation categories.
    pub fn category(self) -> Option<ReasonCategory> {
        match self {
            ConvergedReason::DivergedBreakdown
            | ConvergedReason::DivergedArnoldiRankLoss
            | ConvergedReason::DivergedBreakdownBiCG
            | ConvergedReason::DivergedReducedSystemSingular => Some(ReasonCategory::Breakdown),
            ConvergedReason::DivergedNan => Some(ReasonCategory::Nan),
            ConvergedReason::DivergedInf => Some(ReasonCategory::Inf),
            ConvergedReason::DivergedPcSetupFailed => Some(ReasonCategory::PcSetup),
            ConvergedReason::DivergedPcFailed => Some(ReasonCategory::PcApply),
            _ => None,
        }
    }

    /// PETSc-equivalent convergence reason identifier.
    pub fn petsc_reason(self) -> &'static str {
        match self {
            ConvergedReason::ConvergedRtol => "KSP_CONVERGED_RTOL",
            ConvergedReason::ConvergedAtol => "KSP_CONVERGED_ATOL",
            ConvergedReason::ConvergedTrustRegion => "KSP_CONVERGED_TRUST_REGION",
            ConvergedReason::ConvergedHappyBreakdown => "KSP_CONVERGED_HAPPY_BREAKDOWN",
            ConvergedReason::DivergedNan => "KSP_DIVERGED_NANORINF",
            ConvergedReason::DivergedInf => "KSP_DIVERGED_NANORINF",
            ConvergedReason::DivergedDtol => "KSP_DIVERGED_DTOL",
            ConvergedReason::DivergedMaxIts => "KSP_DIVERGED_ITS",
            ConvergedReason::DivergedBreakdown => "KSP_DIVERGED_BREAKDOWN",
            ConvergedReason::DivergedArnoldiRankLoss => "KSP_DIVERGED_BREAKDOWN",
            ConvergedReason::DivergedBreakdownBiCG => "KSP_DIVERGED_BREAKDOWN_BICG",
            ConvergedReason::DivergedReducedSystemSingular => "KSP_DIVERGED_BREAKDOWN",
            ConvergedReason::DivergedIndefiniteMatrix => "KSP_DIVERGED_INDEFINITE_MAT",
            ConvergedReason::DivergedIndefinitePC => "KSP_DIVERGED_INDEFINITE_PC",
            ConvergedReason::DivergedPcSetupFailed => "KSP_DIVERGED_PCSETUP_FAILED",
            ConvergedReason::DivergedPcFailed => "KSP_DIVERGED_PC_FAILED",
            ConvergedReason::StoppedByMonitor => "KSP_DIVERGED_USER",
            ConvergedReason::Continued => "KSP_CONVERGED_ITERATING",
        }
    }

    /// Classify a non-finite value as a dedicated convergence reason.
    pub fn from_non_finite(value: R) -> Option<Self> {
        if value.is_nan() {
            Some(ConvergedReason::from_failure_kind(FailureReasonKind::Nan))
        } else if value.is_infinite() {
            Some(ConvergedReason::from_failure_kind(FailureReasonKind::Inf))
        } else {
            None
        }
    }

    /// Whether the reason indicates a converged solve.
    pub fn is_converged(self) -> bool {
        matches!(
            self,
            ConvergedReason::ConvergedRtol
                | ConvergedReason::ConvergedAtol
                | ConvergedReason::ConvergedTrustRegion
                | ConvergedReason::ConvergedHappyBreakdown
        )
    }

    /// Whether the reason indicates divergence.
    pub fn is_diverged(self) -> bool {
        !matches!(self, ConvergedReason::Continued) && !self.is_converged()
    }
}

/// Lightweight helper for emitting stable convergence reasons on hot paths.
pub struct ReasonEmitter;

impl ReasonEmitter {
    #[inline(always)]
    pub fn non_finite(value: R) -> Option<ConvergedReason> {
        ConvergedReason::from_non_finite(value)
    }

    #[inline(always)]
    pub fn breakdown() -> ConvergedReason {
        ConvergedReason::from_failure_kind(FailureReasonKind::Breakdown)
    }

    #[inline(always)]
    pub fn breakdown_bicg() -> ConvergedReason {
        ConvergedReason::from_failure_kind(FailureReasonKind::BreakdownBiCG)
    }

    #[inline(always)]
    pub fn indefinite_matrix() -> ConvergedReason {
        ConvergedReason::from_failure_kind(FailureReasonKind::IndefiniteMatrix)
    }

    #[inline(always)]
    pub fn indefinite_pc() -> ConvergedReason {
        ConvergedReason::from_failure_kind(FailureReasonKind::IndefinitePc)
    }

    #[inline]
    pub fn from_error(err: &KError, stage: FailureStage) -> Option<ConvergedReason> {
        map_kerror_to_reason(err, stage)
    }

    #[inline]
    pub fn nested_pc_failure(err: &KError, stage: FailureStage) -> Option<NestedPcFailure> {
        match err {
            KError::PcFailed(msg) => {
                let reason = match stage {
                    FailureStage::Setup => {
                        ConvergedReason::from_failure_kind(FailureReasonKind::PcSetup)
                    }
                    FailureStage::Solve => {
                        ConvergedReason::from_failure_kind(FailureReasonKind::PcApply)
                    }
                };
                Some(NestedPcFailure {
                    component: "pc",
                    reason,
                    iterations: 0,
                    final_norm: None,
                    residual_history_summary: None,
                    detail: format!("stage={stage:?} detail={msg}"),
                })
            }
            KError::FactorError(msg) => {
                let reason = match stage {
                    FailureStage::Setup => {
                        ConvergedReason::from_failure_kind(FailureReasonKind::PcSetup)
                    }
                    FailureStage::Solve => {
                        ConvergedReason::from_failure_kind(FailureReasonKind::PcApply)
                    }
                };
                Some(NestedPcFailure {
                    component: "factorization",
                    reason,
                    iterations: 0,
                    final_norm: None,
                    residual_history_summary: None,
                    detail: format!("stage={stage:?} detail={msg}"),
                })
            }
            KError::ZeroPivot(row) => {
                let reason = match stage {
                    FailureStage::Setup => {
                        ConvergedReason::from_failure_kind(FailureReasonKind::PcSetup)
                    }
                    FailureStage::Solve => {
                        ConvergedReason::from_failure_kind(FailureReasonKind::PcApply)
                    }
                };
                Some(NestedPcFailure {
                    component: "factorization",
                    reason,
                    iterations: 0,
                    final_norm: None,
                    residual_history_summary: None,
                    detail: format!("stage={stage:?} zero_pivot_row={row}"),
                })
            }
            KError::NestedPcFailed(failure) => Some(failure.clone()),
            _ => None,
        }
    }
}

/// Transport-layer mapping of setup/solve errors into stable reason codes.
pub fn map_kerror_to_reason(err: &KError, stage: FailureStage) -> Option<ConvergedReason> {
    match stage {
        FailureStage::Setup => match err {
            KError::NestedPcFailed(failure) => Some(failure.reason),
            KError::PcFailed(_) | KError::FactorError(_) | KError::ZeroPivot(_) => Some(
                ConvergedReason::from_failure_kind(FailureReasonKind::PcSetup),
            ),
            KError::IndefinitePreconditioner | KError::DivergedIndefinitePC => Some(
                ConvergedReason::from_failure_kind(FailureReasonKind::IndefinitePc),
            ),
            KError::BreakdownOrIndefinite => Some(ConvergedReason::from_failure_kind(
                FailureReasonKind::Breakdown,
            )),
            KError::IndefiniteMatrix => Some(ConvergedReason::from_failure_kind(
                FailureReasonKind::IndefiniteMatrix,
            )),
            _ => None,
        },
        FailureStage::Solve => match err {
            KError::NestedPcFailed(failure) => Some(failure.reason),
            KError::PcFailed(_) | KError::FactorError(_) | KError::ZeroPivot(_) => Some(
                ConvergedReason::from_failure_kind(FailureReasonKind::PcApply),
            ),
            KError::BreakdownOrIndefinite => Some(ConvergedReason::from_failure_kind(
                FailureReasonKind::Breakdown,
            )),
            KError::IndefiniteMatrix => Some(ConvergedReason::from_failure_kind(
                FailureReasonKind::IndefiniteMatrix,
            )),
            KError::IndefinitePreconditioner | KError::DivergedIndefinitePC => Some(
                ConvergedReason::from_failure_kind(FailureReasonKind::IndefinitePc),
            ),
            _ => None,
        },
    }
}

/// Counters grouped by stable reason category.
#[derive(Clone, Debug, Default)]
pub struct ReasonDiagnosticsCounters {
    pub breakdown: usize,
    pub nan: usize,
    pub inf: usize,
    pub pc_setup: usize,
    pub pc_apply: usize,
}

impl ReasonDiagnosticsCounters {
    pub fn record_reason(&mut self, reason: ConvergedReason) {
        match reason.category() {
            Some(ReasonCategory::Breakdown) => self.breakdown += 1,
            Some(ReasonCategory::Nan) => self.nan += 1,
            Some(ReasonCategory::Inf) => self.inf += 1,
            Some(ReasonCategory::PcSetup) => self.pc_setup += 1,
            Some(ReasonCategory::PcApply) => self.pc_apply += 1,
            None => {}
        }
    }
}

impl std::fmt::Display for ConvergedReason {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        f.write_str(self.petsc_reason())
    }
}

/// Statistics from a solve operation.
#[derive(Clone, Debug, Default)]
pub struct SolverCounters {
    /// Number of global reduction operations executed by the solver.
    pub num_global_reductions: usize,
    /// Number of reductions that were launched in overlap-friendly/nonblocking paths.
    pub overlap_global_reductions: usize,
    /// Number of residual replacement events performed during the solve.
    pub residual_replacements: usize,
}

/// FGMRES-specific accounting captured by flexible GMRES variants.
#[derive(Clone, Debug, Default)]
pub struct FgmresCounters {
    /// Number of restart boundaries crossed (excluding the initial cycle).
    pub restart_count: usize,
    /// Number of inner Arnoldi iterations completed in the last cycle.
    pub inner_iterations_last_cycle: usize,
    /// Number of orthogonalization passes executed (including second-pass reorthogonalization).
    pub orthog_passes: usize,
    /// Number of happy-breakdown events observed.
    pub happy_breakdowns: usize,
    /// Number of explicit true-residual recomputations (`||b - A x||_2` checks).
    pub explicit_residual_checks: usize,
    /// Number of pipelined-mode fallback/restart interventions triggered by stagnation policy.
    pub pipeline_fallbacks: usize,
    /// Number of mutable preconditioner modification callback invocations.
    pub modify_pc_calls: usize,
    /// Number of pipelined reductions whose completion check was deferred until after local work.
    pub deferred_pipeline_waits: usize,
    /// Number of pipelined reductions that were already complete at finalization time.
    pub immediate_pipeline_completions: usize,
}

/// Estimated reduction counts grouped by algorithm phase.
#[derive(Clone, Debug, Default, PartialEq, Eq)]
pub struct ReductionPhaseDiagnostics {
    /// One-time reductions before entering the main Krylov loop.
    pub startup: usize,
    /// Reductions attributed to the iterative body.
    pub iterative: usize,
    /// Finalization/restart tail reductions after the loop body.
    pub tail: usize,
}

impl SolverCounters {
    /// Infer reduction phase diagnostics from observed totals and an optional model.
    ///
    /// When the model is unavailable, this returns `None`.
    pub fn reduction_phase_diagnostics(
        &self,
        model: Option<&ReductionModel>,
        iterations: usize,
    ) -> Option<ReductionPhaseDiagnostics> {
        let model = model?;
        let startup = model.startup;
        let iterative = (model.per_iteration * iterations as f64).ceil() as usize;
        let modeled_total = startup + iterative + model.tail;
        let tail = if self.num_global_reductions >= startup + iterative {
            self.num_global_reductions - startup - iterative
        } else {
            model.tail
        };
        let _ = modeled_total;
        Some(ReductionPhaseDiagnostics {
            startup,
            iterative,
            tail,
        })
    }
}

/// GCR-specific accounting captured by GCR/PipeGCR variants.
#[derive(Clone, Debug, Default)]
pub struct GcrCounters {
    /// Number of `(p_i, Ap_i)` basis pairs accepted into the Krylov basis.
    pub basis_updates: usize,
    /// Solver-observed synchronization/reduction count.
    pub sync_count: usize,
    /// Number of restart boundaries crossed during the solve.
    pub restart_count: usize,
    /// Whether the solve executed at least one restart.
    pub restarted: bool,
}

/// Structured failure details propagated from a nested preconditioner solve.
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct NestedPcFailure {
    /// Stage/component that produced the nested solve failure.
    pub component: &'static str,
    /// Inner KSP convergence reason.
    pub reason: ConvergedReason,
    /// Inner iteration count.
    pub iterations: usize,
    /// Final inner norm semantics/summary reported by the nested solve.
    pub final_norm: Option<String>,
    /// Inner residual history summary captured by nested monitors or solver callbacks.
    pub residual_history_summary: Option<String>,
    /// Human-readable context (inner solver/pc selections, etc.).
    pub detail: String,
}

#[cfg(feature = "metrics")]
#[derive(Clone, Debug, Default)]
pub struct SolveMetrics {
    pub reductions: usize,
    pub reduction_wait_nanos: u64,
    pub reduction_overlap_nanos: u64,
    pub matvec_nanos: u64,
    pub pc_apply_nanos: u64,
    pub bytes_reduced: usize,
}

#[cfg(not(feature = "metrics"))]
#[derive(Clone, Debug, Default)]
pub struct SolveMetrics;

/// Reduction accounting model for a solver/variant.
#[derive(Clone, Debug)]
pub struct ReductionModel {
    /// Variant/algorithm label used in diagnostics output.
    pub variant: &'static str,
    /// One-time reductions before the main iteration loop.
    pub startup: usize,
    /// Typical global reductions performed per iteration.
    pub per_iteration: f64,
    /// Additional reductions paid on finalize/restart paths.
    pub tail: usize,
}

impl ReductionModel {
    /// Conservative estimate for `num_global_reductions` given observed iterations.
    pub fn estimate_total(&self, iterations: usize) -> usize {
        self.startup + (self.per_iteration * iterations as f64).ceil() as usize + self.tail
    }
}

/// Statistics from a solve operation.
#[must_use]
#[derive(Clone, Debug)]
pub struct SolveStats<R> {
    /// Number of iterations performed
    pub iterations: usize,
    /// Canonical final residual norm for reporting.
    ///
    /// Solvers should prefer storing the explicit true residual norm when it is
    /// computed; otherwise this may be a recurrence/estimated residual.
    pub final_residual: R,
    /// Optional final residual from the Krylov recurrence/estimator.
    pub final_recurrence_residual: Option<R>,
    /// Optional final explicit true residual `||b - A x||_2`.
    pub final_true_residual: Option<R>,
    /// Optional latest preconditioned residual norm.
    pub last_preconditioned_residual: Option<R>,
    /// Reason for stopping.
    ///
    /// This includes mapped preconditioner lifecycle failures (for example
    /// `DivergedPcSetupFailed`/`DivergedPcFailed`) when setup/apply errors are
    /// propagated through `KspContext` into `SolveStats`.
    pub reason: ConvergedReason,
    /// Additional counters collected during the solve.
    pub counters: SolverCounters,
    /// Optional GCR-specific counters.
    pub gcr_counters: Option<GcrCounters>,
    /// Optional FGMRES-specific counters.
    pub fgmres_counters: Option<FgmresCounters>,
    /// Reduction accounting model for the selected solver variant.
    pub reduction_model: Option<ReductionModel>,
    /// Total number of complex drift events observed during reductions.
    pub complex_drift_events: usize,
    /// Per-kind complex drift counts captured by the solver.
    pub complex_drift_counts: [usize; 6],
    /// Maximum relative imaginary magnitude observed.
    pub complex_drift_max_rel: R,
    /// Optional solver timing and reduction metrics.
    pub metrics: SolveMetrics,
    /// Structured details for nested preconditioner failures, when available.
    ///
    /// Outer solves use this to preserve the inner component/reason/iteration
    /// context when KSP-as-PC or other nested preconditioners fail.
    pub nested_pc_failure: Option<NestedPcFailure>,
    /// Stable per-category reason counters for diagnostics automation.
    pub reason_counters: ReasonDiagnosticsCounters,
    /// True when GMRES internally retried by resetting `x` and falling back to
    /// a classical GMRES pass.
    pub gmres_classical_retry: bool,
    /// Contract acceptance status from true residual + stop reason classification.
    pub acceptance_status: AcceptanceStatus,
    /// Breakdown/divergence reason captured before any residual-based override.
    pub breakdown_reason: Option<ConvergedReason>,
    /// Optional note describing residual-based acceptance override decisions.
    pub residual_override_note: Option<String>,
    /// Count of orthogonalization passes executed (including reorthogonalization).
    pub orthogonalization_passes: usize,
    /// Whether an orthogonalization rank-loss failure was observed.
    pub orthogonalization_rank_loss: bool,
    /// Max estimate of orthogonality loss observed during the solve.
    pub max_orthogonality_loss_estimate: R,
    /// Effective solver variant selected at runtime (if policy hooks mutated it).
    pub effective_variant: Option<String>,
    /// Effective restart value selected at runtime.
    pub effective_restart: Option<usize>,
    /// Effective residual-check policy selected at runtime.
    pub effective_residual_check_policy: Option<String>,
}

impl<R: Default> SolveStats<R> {
    /// Construct a new statistics record with zeroed counters.
    pub fn new(iterations: usize, final_residual: R, reason: ConvergedReason) -> Self {
        Self {
            iterations,
            final_residual,
            final_recurrence_residual: None,
            final_true_residual: None,
            last_preconditioned_residual: None,
            reason,
            counters: SolverCounters::default(),
            gcr_counters: None,
            fgmres_counters: None,
            reduction_model: None,
            complex_drift_events: 0,
            complex_drift_counts: [0; 6],
            complex_drift_max_rel: R::default(),
            metrics: SolveMetrics::default(),
            nested_pc_failure: None,
            reason_counters: ReasonDiagnosticsCounters::default(),
            gmres_classical_retry: false,
            acceptance_status: if reason.is_converged() {
                AcceptanceStatus::Ok
            } else {
                AcceptanceStatus::Failed
            },
            breakdown_reason: None,
            residual_override_note: None,
            orthogonalization_passes: 0,
            orthogonalization_rank_loss: false,
            max_orthogonality_loss_estimate: R::default(),
            effective_variant: None,
            effective_restart: None,
            effective_residual_check_policy: None,
        }
    }

    /// Attach solver counters to an existing statistics record.
    pub fn with_counters(mut self, counters: SolverCounters) -> Self {
        self.counters = counters;
        self
    }

    pub fn with_reduction_model(mut self, model: ReductionModel) -> Self {
        self.reduction_model = Some(model);
        self
    }

    /// Attach GCR-specific counters.
    pub fn with_gcr_counters(mut self, counters: GcrCounters) -> Self {
        self.gcr_counters = Some(counters);
        self
    }

    /// Attach FGMRES-specific counters.
    pub fn with_fgmres_counters(mut self, counters: FgmresCounters) -> Self {
        self.fgmres_counters = Some(counters);
        self
    }

    /// Attach nested preconditioner failure metadata.
    pub fn with_nested_pc_failure(mut self, failure: NestedPcFailure) -> Self {
        self.nested_pc_failure = Some(failure);
        self
    }

    /// Mark whether GMRES performed an internal retry to classical GMRES.
    pub fn with_gmres_classical_retry(mut self, used: bool) -> Self {
        self.gmres_classical_retry = used;
        self
    }

    /// Finalize stable reason diagnostics for this stats record.
    pub fn finalize_reason_counters(mut self) -> Self {
        self.reason_counters.record_reason(self.reason);
        if let Some(inner) = self.nested_pc_failure.as_ref() {
            self.reason_counters.record_reason(inner.reason);
        }
        self
    }

    /// Per-phase reduction diagnostics inferred from solver counters/model.
    pub fn reduction_phase_diagnostics(&self) -> Option<ReductionPhaseDiagnostics> {
        self.counters
            .reduction_phase_diagnostics(self.reduction_model.as_ref(), self.iterations)
    }

    /// Attach orthogonalization diagnostics.
    pub fn with_orthogonalization_diagnostics(
        mut self,
        passes: usize,
        rank_loss: bool,
        max_loss_estimate: R,
    ) -> Self {
        self.orthogonalization_passes = passes;
        self.orthogonalization_rank_loss = rank_loss;
        self.max_orthogonality_loss_estimate = max_loss_estimate;
        self
    }

    /// Attach runtime-effective policy knobs selected by the solver.
    pub fn with_effective_runtime_policy(
        mut self,
        variant: impl Into<String>,
        restart: usize,
        residual_check_policy: impl Into<String>,
    ) -> Self {
        self.effective_variant = Some(variant.into());
        self.effective_restart = Some(restart);
        self.effective_residual_check_policy = Some(residual_check_policy.into());
        self
    }
}

impl Convergence {
    /// Create new convergence criteria.
    pub fn new(rtol: R, atol: R, dtol: R, max_iters: usize) -> Self {
        Self {
            rtol,
            atol,
            dtol,
            max_iters,
        }
    }

    /// Check convergence/divergence criteria.
    ///
    /// Returns (reason, SolveStats) based on current residual norm and iteration count.
    ///
    /// # Arguments
    /// * `rnorm` - Current residual norm ‖r‖
    /// * `bnorm` - Right-hand side norm ‖b‖
    /// * `iters` - Current iteration count
    ///
    /// # Returns
    /// Tuple of (ConvergedReason, SolveStats) indicating the stopping reason.
    pub fn check(&self, rnorm: R, bnorm: R, iters: usize) -> (ConvergedReason, SolveStats<R>) {
        if let Some(reason) = ConvergedReason::from_non_finite(rnorm)
            .or_else(|| ConvergedReason::from_non_finite(bnorm))
        {
            let stats = SolveStats::new(iters, rnorm, reason);
            return (reason, stats);
        }

        // Absolute tolerance test first (most restrictive)
        if rnorm <= self.atol {
            let stats = SolveStats::new(iters, rnorm, ConvergedReason::ConvergedAtol);
            return (ConvergedReason::ConvergedAtol, stats);
        }

        // Relative tolerance test
        if rnorm <= self.rtol * bnorm {
            let stats = SolveStats::new(iters, rnorm, ConvergedReason::ConvergedRtol);
            return (ConvergedReason::ConvergedRtol, stats);
        }

        // Divergence test
        if rnorm >= self.dtol * bnorm {
            let stats = SolveStats::new(iters, rnorm, ConvergedReason::DivergedDtol);
            return (ConvergedReason::DivergedDtol, stats);
        }

        // Maximum iterations test
        if iters >= self.max_iters {
            let stats = SolveStats::new(iters, rnorm, ConvergedReason::DivergedMaxIts);
            return (ConvergedReason::DivergedMaxIts, stats);
        }

        // Continue iterating
        let stats = SolveStats::new(iters, rnorm, ConvergedReason::Continued);
        (ConvergedReason::Continued, stats)
    }
}

// Legacy convenience method for backward compatibility
impl Convergence {
    /// Legacy method for backward compatibility.
    /// Returns (should_stop, stats) given current `res_norm` and iteration `i`.
    ///
    /// **Deprecated**: Use `check()` instead for more detailed convergence information.
    #[deprecated(since = "0.1.0", note = "use check() method instead")]
    pub fn check_legacy(&self, res_norm: R, res0_norm: R, i: usize) -> (bool, SolveStats<R>) {
        let (reason, stats) = self.check(res_norm, res0_norm, i);
        let converged = matches!(
            reason,
            ConvergedReason::ConvergedRtol | ConvergedReason::ConvergedAtol
        );
        let mut legacy_stats =
            SolveStats::new(stats.iterations, stats.final_residual, stats.reason);
        legacy_stats.counters = stats.counters;
        (
            converged || reason != ConvergedReason::Continued,
            legacy_stats,
        )
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_convergence_new() {
        let conv = Convergence::new(1e-6, 1e-12, 1e3, 1000);
        assert_eq!(conv.rtol, 1e-6);
        assert_eq!(conv.atol, 1e-12);
        assert_eq!(conv.dtol, 1e3);
        assert_eq!(conv.max_iters, 1000);
    }

    #[test]
    fn test_converged_absolute_tolerance() {
        let conv = Convergence::new(1e-6, 1e-8, 1e3, 100);
        let rnorm = 1e-9; // Less than atol
        let bnorm = 1.0;
        let iters = 5;

        let (reason, stats) = conv.check(rnorm, bnorm, iters);

        assert_eq!(reason, ConvergedReason::ConvergedAtol);
        assert_eq!(stats.reason, ConvergedReason::ConvergedAtol);
        assert_eq!(stats.iterations, 5);
        assert_eq!(stats.final_residual, 1e-9);
    }

    #[test]
    fn test_converged_relative_tolerance() {
        let conv = Convergence::new(1e-6, 1e-12, 1e3, 100);
        let rnorm = 1e-7; // Greater than atol but satisfies rtol
        let bnorm = 1.0;
        let iters = 10;

        let (reason, stats) = conv.check(rnorm, bnorm, iters);

        assert_eq!(reason, ConvergedReason::ConvergedRtol);
        assert_eq!(stats.reason, ConvergedReason::ConvergedRtol);
        assert_eq!(stats.iterations, 10);
        assert_eq!(stats.final_residual, 1e-7);
    }

    #[test]
    fn test_diverged_tolerance() {
        let conv = Convergence::new(1e-6, 1e-12, 2.0, 100);
        let rnorm = 3.0; // Greater than dtol * bnorm
        let bnorm = 1.0;
        let iters = 5;

        let (reason, stats) = conv.check(rnorm, bnorm, iters);

        assert_eq!(reason, ConvergedReason::DivergedDtol);
        assert_eq!(stats.reason, ConvergedReason::DivergedDtol);
        assert_eq!(stats.iterations, 5);
        assert_eq!(stats.final_residual, 3.0);
    }

    #[test]
    fn test_diverged_max_iterations() {
        let conv = Convergence::new(1e-6, 1e-12, 1e3, 10);
        let rnorm = 1e-3; // Not converged but within tolerance bounds
        let bnorm = 1.0;
        let iters = 10; // Equal to max_iters

        let (reason, stats) = conv.check(rnorm, bnorm, iters);

        assert_eq!(reason, ConvergedReason::DivergedMaxIts);
        assert_eq!(stats.reason, ConvergedReason::DivergedMaxIts);
        assert_eq!(stats.iterations, 10);
        assert_eq!(stats.final_residual, 1e-3);
    }

    #[test]
    fn test_continued() {
        let conv = Convergence::new(1e-6, 1e-12, 1e3, 100);
        let rnorm = 1e-3; // Not converged, not diverged, within iteration limit
        let bnorm = 1.0;
        let iters = 5;

        let (reason, stats) = conv.check(rnorm, bnorm, iters);

        assert_eq!(reason, ConvergedReason::Continued);
        assert_eq!(stats.reason, ConvergedReason::Continued);
        assert_eq!(stats.iterations, 5);
        assert_eq!(stats.final_residual, 1e-3);
    }

    #[test]
    fn test_convergence_precedence() {
        // Absolute tolerance takes precedence over relative
        let conv = Convergence::new(1e-6, 1e-8, 1e3, 100);
        let rnorm = 1e-9; // Satisfies both atol and rtol
        let bnorm = 1.0;
        let iters = 5;

        let (reason, _) = conv.check(rnorm, bnorm, iters);
        assert_eq!(reason, ConvergedReason::ConvergedAtol);
    }

    #[test]
    fn test_converged_reason_equality() {
        assert_eq!(
            ConvergedReason::ConvergedRtol,
            ConvergedReason::ConvergedRtol
        );
        assert_eq!(
            ConvergedReason::ConvergedAtol,
            ConvergedReason::ConvergedAtol
        );
        assert_eq!(ConvergedReason::DivergedDtol, ConvergedReason::DivergedDtol);
        assert_eq!(
            ConvergedReason::DivergedMaxIts,
            ConvergedReason::DivergedMaxIts
        );
        assert_eq!(ConvergedReason::Continued, ConvergedReason::Continued);

        assert_ne!(
            ConvergedReason::ConvergedRtol,
            ConvergedReason::ConvergedAtol
        );
        assert_ne!(
            ConvergedReason::DivergedDtol,
            ConvergedReason::DivergedMaxIts
        );
    }

    #[test]
    fn test_converged_reason_debug() {
        let reason = ConvergedReason::ConvergedRtol;
        let debug_str = format!("{:?}", reason);
        assert!(debug_str.contains("ConvergedRtol"));
    }

    #[test]
    fn test_convergence_nan_or_inf() {
        let conv = Convergence::new(1e-6, 1e-12, 1e3, 10);
        let (reason_nan, _) = conv.check(f64::NAN, 1.0, 1);
        assert_eq!(reason_nan, ConvergedReason::DivergedNan);

        let (reason_inf, _) = conv.check(f64::INFINITY, 1.0, 2);
        assert_eq!(reason_inf, ConvergedReason::DivergedInf);

        let (reason_bnorm_inf, _) = conv.check(1.0, f64::INFINITY, 3);
        assert_eq!(reason_bnorm_inf, ConvergedReason::DivergedInf);
    }

    #[test]
    fn test_solve_stats_clone() {
        let stats = SolveStats::new(42, 1e-8, ConvergedReason::ConvergedRtol);

        let cloned = stats.clone();
        assert_eq!(cloned.iterations, 42);
        assert_eq!(cloned.final_residual, 1e-8);
        assert_eq!(cloned.final_recurrence_residual, None);
        assert_eq!(cloned.final_true_residual, None);
        assert_eq!(cloned.last_preconditioned_residual, None);
        assert_eq!(cloned.reason, ConvergedReason::ConvergedRtol);
    }

    #[test]
    fn test_solve_stats_debug() {
        let stats = SolveStats::new(10, 1e-6, ConvergedReason::ConvergedAtol);

        let debug_str = format!("{:?}", stats);
        assert!(debug_str.contains("10"));
        assert!(debug_str.contains("ConvergedAtol"));
    }

    #[test]
    #[allow(deprecated)]
    fn test_legacy_check_convergence() {
        let conv = Convergence::new(1e-6, 1e-12, 1e3, 100);
        let res_norm = 1e-8;
        let res0_norm = 1.0;
        let iters = 5;

        let (should_stop, stats) = conv.check_legacy(res_norm, res0_norm, iters);

        assert!(should_stop);
        assert_eq!(stats.iterations, 5);
        assert_eq!(stats.final_residual, 1e-8);
    }

    #[test]
    #[allow(deprecated)]
    fn test_legacy_check_continue() {
        let conv = Convergence::new(1e-6, 1e-12, 1e3, 100);
        let res_norm = 1e-3;
        let res0_norm = 1.0;
        let iters = 5;

        let (should_stop, stats) = conv.check_legacy(res_norm, res0_norm, iters);

        assert!(!should_stop);
        assert_eq!(stats.iterations, 5);
        assert_eq!(stats.final_residual, 1e-3);
        assert_eq!(stats.reason, ConvergedReason::Continued);
    }

    #[test]
    fn nested_pc_failure_mapping_preserves_structured_inner_reason() {
        let err = KError::NestedPcFailed(NestedPcFailure {
            component: "pc_ksp",
            reason: ConvergedReason::DivergedBreakdown,
            iterations: 4,
            final_norm: Some("true_residual_l2=1.0e+00".into()),
            residual_history_summary: Some("history_len=4".into()),
            detail: "inner failure".into(),
        });
        let reason = map_kerror_to_reason(&err, FailureStage::Solve).expect("reason");
        assert_eq!(reason, ConvergedReason::DivergedBreakdown);
        let nested =
            ReasonEmitter::nested_pc_failure(&err, FailureStage::Solve).expect("nested metadata");
        assert_eq!(nested.component, "pc_ksp");
        assert_eq!(nested.reason, ConvergedReason::DivergedBreakdown);
        assert_eq!(nested.iterations, 4);
    }

    #[test]
    fn test_different_numeric_types() {
        // Test with f64 (which implements From<f64>)
        let conv_f64 = Convergence::new(1e-6f64, 1e-12f64, 1e3f64, 100);
        let (reason, _) = conv_f64.check(1e-8f64, 1.0f64, 5);
        assert_eq!(reason, ConvergedReason::ConvergedRtol);

        // Test with different tolerances
        let conv2 = Convergence::new(1e-8, 1e-16, 1e6, 50);
        let (reason2, _) = conv2.check(1e-10, 1.0, 10);
        assert_eq!(reason2, ConvergedReason::ConvergedRtol);
    }
}