ruvector_solver/

neumann.rs

//! Jacobi-preconditioned Neumann Series iterative solver.
//!
//! Solves the linear system `Ax = b` by splitting `A = D - R` (where `D` is
//! the diagonal part) and iterating:
//!
//! ```text
//! x_{k+1} = x_k + D^{-1} (b - A x_k)
//! ```
//!
//! This is equivalent to the Neumann series `x = sum_{k=0}^{K} M^k D^{-1} b`
//! where `M = I - D^{-1} A`. Convergence requires `rho(M) < 1`, which is
//! guaranteed for strictly diagonally dominant matrices.
//!
//! # Algorithm
//!
//! The iteration maintains a running solution `x` and residual `r = b - Ax`:
//!
//! ```text
//! x_0 = D^{-1} b
//! for k = 0, 1, 2, ...:
//!     r = b - A * x_k
//!     x_{k+1} = x_k + D^{-1} * r
//!     if ||r|| < tolerance:
//!         break
//! ```
//!
//! # Convergence
//!
//! Before solving, the solver estimates `rho(I - D^{-1}A)` via a 10-step
//! power iteration and rejects the problem with
//! [`SolverError::SpectralRadiusExceeded`] if `rho >= 1.0`. During iteration,
//! if the residual grows by more than 2x between consecutive steps,
//! [`SolverError::NumericalInstability`] is returned.

use std::time::Instant;

use tracing::{debug, info, instrument, warn};

use crate::error::{SolverError, ValidationError};
use crate::traits::SolverEngine;
use crate::types::{
    Algorithm, ComplexityClass, ComplexityEstimate, ComputeBudget, ConvergenceInfo, CsrMatrix,
    SolverResult, SparsityProfile,
};

/// Number of power-iteration steps used to estimate the spectral radius.
const POWER_ITERATION_STEPS: usize = 10;

/// If the residual grows by more than this factor in a single step, the solver
/// declares numerical instability.
const INSTABILITY_GROWTH_FACTOR: f64 = 2.0;

// ---------------------------------------------------------------------------
// NeumannSolver
// ---------------------------------------------------------------------------

/// Neumann Series solver for sparse linear systems.
///
/// Computes `x = sum_{k=0}^{K} (I - A)^k * b` by maintaining a residual
/// vector and accumulating partial sums until convergence.
///
/// # Example
///
/// ```rust
/// use ruvector_solver::types::CsrMatrix;
/// use ruvector_solver::neumann::NeumannSolver;
///
/// // Diagonally dominant 2x2: A = [[2, -0.5], [-0.5, 2]]
/// let a = CsrMatrix::<f32>::from_coo(2, 2, vec![
///     (0, 0, 2.0_f32), (0, 1, -0.5_f32),
///     (1, 0, -0.5_f32), (1, 1, 2.0_f32),
/// ]);
/// let b = vec![1.0_f32, 1.0];
///
/// let solver = NeumannSolver::new(1e-6, 500);
/// let result = solver.solve(&a, &b).unwrap();
/// assert!(result.residual_norm < 1e-4);
/// ```
#[derive(Debug, Clone)]
pub struct NeumannSolver {
    /// Target residual L2 norm for convergence.
    pub tolerance: f64,
    /// Maximum number of iterations before giving up.
    pub max_iterations: usize,
}

impl NeumannSolver {
    /// Create a new `NeumannSolver`.
    ///
    /// # Arguments
    ///
    /// * `tolerance` - Stop when `||r|| < tolerance`.
    /// * `max_iterations` - Upper bound on iterations.
    pub fn new(tolerance: f64, max_iterations: usize) -> Self {
        Self {
            tolerance,
            max_iterations,
        }
    }

    /// Estimate the spectral radius of `M = I - D^{-1}A` via 10-step power
    /// iteration.
    ///
    /// Runs [`POWER_ITERATION_STEPS`] iterations of the power method on the
    /// Jacobi iteration matrix `M = I - D^{-1}A`. Returns the Rayleigh-quotient
    /// estimate of the dominant eigenvalue magnitude.
    ///
    /// # Arguments
    ///
    /// * `matrix` - The coefficient matrix `A` (must be square).
    ///
    /// # Returns
    ///
    /// Estimated `|lambda_max(I - D^{-1}A)|`. If this is `>= 1.0`, the
    /// Jacobi-preconditioned Neumann series will diverge.
    #[instrument(skip(matrix), fields(n = matrix.rows))]
    pub fn estimate_spectral_radius(matrix: &CsrMatrix<f32>) -> f64 {
        let n = matrix.rows;
        if n == 0 {
            return 0.0;
        }

        let d_inv = extract_diag_inv_f32(matrix);
        Self::estimate_spectral_radius_with_diag(matrix, &d_inv)
    }

    /// Inner helper: estimate spectral radius using a pre-computed `d_inv`.
    ///
    /// This avoids recomputing the diagonal inverse when the caller already
    /// has it (e.g. `solve()` needs `d_inv` for both the spectral check and
    /// the Jacobi iteration).
    fn estimate_spectral_radius_with_diag(
        matrix: &CsrMatrix<f32>,
        d_inv: &[f32],
    ) -> f64 {
        let n = matrix.rows;
        if n == 0 {
            return 0.0;
        }

        // Initialise with a deterministic pseudo-random unit vector.
        let mut v: Vec<f32> = (0..n)
            .map(|i| ((i.wrapping_mul(7).wrapping_add(13)) % 100) as f32 / 100.0)
            .collect();
        let norm = l2_norm_f32(&v);
        if norm > 1e-12 {
            scale_vec_f32(&mut v, 1.0 / norm);
        }

        let mut av = vec![0.0f32; n]; // scratch for A*v
        let mut w = vec![0.0f32; n]; // scratch for M*v = v - D^{-1}*A*v
        let mut eigenvalue_estimate = 0.0f64;

        for _ in 0..POWER_ITERATION_STEPS {
            // w = v - D^{-1} * A * v  (i.e. M * v)
            matrix.spmv(&v, &mut av);
            for j in 0..n {
                w[j] = v[j] - d_inv[j] * av[j];
            }

            // Rayleigh quotient: lambda = v^T w  (v is unit-length).
            let dot: f64 = v
                .iter()
                .zip(w.iter())
                .map(|(&a, &b)| a as f64 * b as f64)
                .sum();
            eigenvalue_estimate = dot;

            // Normalise w -> v for the next step.
            let w_norm = l2_norm_f32(&w);
            if w_norm < 1e-12 {
                break;
            }
            for j in 0..n {
                v[j] = w[j] / w_norm as f32;
            }
        }

        let rho = eigenvalue_estimate.abs();
        debug!(rho, "estimated spectral radius of (I - D^-1 A)");
        rho
    }

    /// Core Jacobi-preconditioned Neumann-series solve operating on `f32`.
    ///
    /// Validates inputs, checks the spectral radius of `I - D^{-1}A` via
    /// power iteration, then runs the iteration returning a [`SolverResult`].
    ///
    /// # Errors
    ///
    /// - [`SolverError::InvalidInput`] if the matrix is non-square or the RHS
    ///   length does not match.
    /// - [`SolverError::SpectralRadiusExceeded`] if `rho(I - D^{-1}A) >= 1`.
    /// - [`SolverError::NumericalInstability`] if the residual grows by more
    ///   than 2x in a single step.
    /// - [`SolverError::NonConvergence`] if the iteration budget is exhausted.
    #[instrument(skip(self, matrix, rhs), fields(n = matrix.rows, nnz = matrix.nnz()))]
    pub fn solve(
        &self,
        matrix: &CsrMatrix<f32>,
        rhs: &[f32],
    ) -> Result<SolverResult, SolverError> {
        let start = Instant::now();
        let n = matrix.rows;

        // ------------------------------------------------------------------
        // Input validation
        // ------------------------------------------------------------------
        if matrix.rows != matrix.cols {
            return Err(SolverError::InvalidInput(
                ValidationError::DimensionMismatch(format!(
                    "matrix must be square: got {}x{}",
                    matrix.rows, matrix.cols,
                )),
            ));
        }

        if rhs.len() != n {
            return Err(SolverError::InvalidInput(
                ValidationError::DimensionMismatch(format!(
                    "rhs length {} does not match matrix dimension {}",
                    rhs.len(),
                    n,
                )),
            ));
        }

        // Edge case: empty system.
        if n == 0 {
            return Ok(SolverResult {
                solution: Vec::new(),
                iterations: 0,
                residual_norm: 0.0,
                wall_time: start.elapsed(),
                convergence_history: Vec::new(),
                algorithm: Algorithm::Neumann,
            });
        }

        // Extract D^{-1} once — reused for both the spectral radius check
        // and the Jacobi-preconditioned iteration that follows.
        let d_inv = extract_diag_inv_f32(matrix);

        // ------------------------------------------------------------------
        // Spectral radius pre-check (10-step power iteration on I - D^{-1}A)
        // ------------------------------------------------------------------
        let rho = Self::estimate_spectral_radius_with_diag(matrix, &d_inv);
        if rho >= 1.0 {
            warn!(rho, "spectral radius >= 1.0, Neumann series will diverge");
            return Err(SolverError::SpectralRadiusExceeded {
                spectral_radius: rho,
                limit: 1.0,
                algorithm: Algorithm::Neumann,
            });
        }
        info!(rho, "spectral radius check passed");

        // ------------------------------------------------------------------
        // Jacobi-preconditioned iteration (fused kernel)
        //
        //   x_0 = D^{-1} * b
        //   loop:
        //       r = b - A * x_k            (fused with norm computation)
        //       if ||r|| < tolerance: break
        //       x_{k+1} = x_k + D^{-1} * r (fused with residual storage)
        //
        // Key optimization: uses fused_residual_norm_sq to compute
        // r = b - Ax and ||r||^2 in a single pass, avoiding a separate
        // spmv + subtraction + norm computation (3 memory traversals -> 1).
        // ------------------------------------------------------------------
        let mut x: Vec<f32> = (0..n).map(|i| d_inv[i] * rhs[i]).collect();
        let mut r = vec![0.0f32; n]; // residual buffer (reused each iteration)

        let mut convergence_history = Vec::with_capacity(self.max_iterations.min(256));
        let mut prev_residual_norm = f64::MAX;
        let final_residual_norm: f64;
        let mut iterations_done: usize = 0;

        for k in 0..self.max_iterations {
            // Fused: compute r = b - Ax and ||r||^2 in one pass.
            let residual_norm_sq = matrix.fused_residual_norm_sq(&x, rhs, &mut r);
            let residual_norm = residual_norm_sq.sqrt();
            iterations_done = k + 1;

            convergence_history.push(ConvergenceInfo {
                iteration: k,
                residual_norm,
            });

            debug!(iteration = k, residual_norm, "neumann iteration");

            // Convergence check.
            if residual_norm < self.tolerance {
                final_residual_norm = residual_norm;
                info!(iterations = iterations_done, residual_norm, "converged");
                return Ok(SolverResult {
                    solution: x,
                    iterations: iterations_done,
                    residual_norm: final_residual_norm,
                    wall_time: start.elapsed(),
                    convergence_history,
                    algorithm: Algorithm::Neumann,
                });
            }

            // NaN / Inf guard.
            if residual_norm.is_nan() || residual_norm.is_infinite() {
                return Err(SolverError::NumericalInstability {
                    iteration: k,
                    detail: format!("residual became {residual_norm}"),
                });
            }

            // Instability check: residual grew by > 2x.
            if k > 0
                && prev_residual_norm < f64::MAX
                && prev_residual_norm > 0.0
                && residual_norm > INSTABILITY_GROWTH_FACTOR * prev_residual_norm
            {
                warn!(
                    iteration = k,
                    prev = prev_residual_norm,
                    current = residual_norm,
                    "residual diverging",
                );
                return Err(SolverError::NumericalInstability {
                    iteration: k,
                    detail: format!(
                        "residual grew from {prev_residual_norm:.6e} to \
                         {residual_norm:.6e} (>{INSTABILITY_GROWTH_FACTOR:.0}x)",
                    ),
                });
            }

            // Fused update: x[j] += d_inv[j] * r[j]
            // 4-wide unrolled for ILP.
            let chunks = n / 4;
            for c in 0..chunks {
                let j = c * 4;
                x[j] += d_inv[j] * r[j];
                x[j + 1] += d_inv[j + 1] * r[j + 1];
                x[j + 2] += d_inv[j + 2] * r[j + 2];
                x[j + 3] += d_inv[j + 3] * r[j + 3];
            }
            for j in (chunks * 4)..n {
                x[j] += d_inv[j] * r[j];
            }

            prev_residual_norm = residual_norm;
        }

        // Exhausted iteration budget without converging.
        final_residual_norm = prev_residual_norm;
        Err(SolverError::NonConvergence {
            iterations: iterations_done,
            residual: final_residual_norm,
            tolerance: self.tolerance,
        })
    }
}

// ---------------------------------------------------------------------------
// SolverEngine trait implementation (f64 interface)
// ---------------------------------------------------------------------------

impl SolverEngine for NeumannSolver {
    /// Solve via the Neumann series.
    ///
    /// Adapts the `f64` trait interface to the internal `f32` solver by
    /// converting the input matrix and RHS, running the solver, then
    /// returning the `f32` solution.
    fn solve(
        &self,
        matrix: &CsrMatrix<f64>,
        rhs: &[f64],
        budget: &ComputeBudget,
    ) -> Result<SolverResult, SolverError> {
        let start = Instant::now();

        // Validate that f64 values fit in f32 range.
        for (i, &v) in matrix.values.iter().enumerate() {
            if v.is_finite() && v.abs() > f32::MAX as f64 {
                return Err(SolverError::InvalidInput(
                    ValidationError::NonFiniteValue(format!(
                        "matrix value at index {i} ({v:.6e}) overflows f32"
                    )),
                ));
            }
        }
        for (i, &v) in rhs.iter().enumerate() {
            if v.is_finite() && v.abs() > f32::MAX as f64 {
                return Err(SolverError::InvalidInput(
                    ValidationError::NonFiniteValue(format!(
                        "rhs value at index {i} ({v:.6e}) overflows f32"
                    )),
                ));
            }
        }

        // Convert f64 matrix to f32 for the core solver.
        let f32_matrix = CsrMatrix {
            row_ptr: matrix.row_ptr.clone(),
            col_indices: matrix.col_indices.clone(),
            values: matrix.values.iter().map(|&v| v as f32).collect(),
            rows: matrix.rows,
            cols: matrix.cols,
        };
        let f32_rhs: Vec<f32> = rhs.iter().map(|&v| v as f32).collect();

        // Use the tighter of the solver's own tolerance and the caller's budget,
        // but no tighter than f32 precision allows (the Neumann solver operates
        // internally in f32, so residuals below ~f32::EPSILON are unreachable).
        let max_iters = self.max_iterations.min(budget.max_iterations);
        let tol = self
            .tolerance
            .min(budget.tolerance)
            .max(f32::EPSILON as f64 * 4.0);

        let inner_solver = NeumannSolver::new(tol, max_iters);

        let mut result = inner_solver.solve(&f32_matrix, &f32_rhs)?;

        // Check wall-time budget.
        if start.elapsed() > budget.max_time {
            return Err(SolverError::BudgetExhausted {
                reason: "wall-clock time limit exceeded".to_string(),
                elapsed: start.elapsed(),
            });
        }

        // Adjust wall time to include conversion overhead.
        result.wall_time = start.elapsed();
        Ok(result)
    }

    fn estimate_complexity(
        &self,
        profile: &SparsityProfile,
        n: usize,
    ) -> ComplexityEstimate {
        // Estimated iterations: ceil( ln(1/tol) / |ln(rho)| )
        let rho = profile.estimated_spectral_radius.max(0.01).min(0.999);
        let est_iters = ((1.0 / self.tolerance).ln() / (1.0 - rho).ln().abs())
            .ceil() as usize;
        let est_iters = est_iters.min(self.max_iterations).max(1);

        ComplexityEstimate {
            algorithm: Algorithm::Neumann,
            // Each iteration does one SpMV (2 * nnz flops) + O(n) vector ops.
            estimated_flops: (est_iters as u64) * (profile.nnz as u64) * 2,
            estimated_iterations: est_iters,
            // Working memory: x, r, ar (3 vectors of f32).
            estimated_memory_bytes: n * 4 * 3,
            complexity_class: ComplexityClass::SublinearNnz,
        }
    }

    fn algorithm(&self) -> Algorithm {
        Algorithm::Neumann
    }
}

// ---------------------------------------------------------------------------
// Internal helpers
// ---------------------------------------------------------------------------

/// Extract `D^{-1}` from a CSR matrix (the reciprocal of each diagonal entry).
///
/// If a diagonal entry is zero or very small, uses `1.0` as a fallback to
/// avoid division by zero.
fn extract_diag_inv_f32(matrix: &CsrMatrix<f32>) -> Vec<f32> {
    let n = matrix.rows;
    let mut d_inv = vec![1.0f32; n];
    for i in 0..n {
        let start = matrix.row_ptr[i];
        let end = matrix.row_ptr[i + 1];
        for idx in start..end {
            if matrix.col_indices[idx] == i {
                let diag = matrix.values[idx];
                if diag.abs() > 1e-15 {
                    d_inv[i] = 1.0 / diag;
                } else {
                    warn!(
                        row = i,
                        diag_value = %diag,
                        "zero or near-zero diagonal entry; substituting 1.0 — matrix may be singular"
                    );
                }
                break;
            }
        }
    }
    d_inv
}

/// Compute the L2 (Euclidean) norm of a slice of `f32` values.
///
/// Uses `f64` accumulation to reduce catastrophic cancellation on large
/// vectors.
#[inline]
fn l2_norm_f32(v: &[f32]) -> f32 {
    let sum: f64 = v.iter().map(|&x| (x as f64) * (x as f64)).sum();
    sum.sqrt() as f32
}

/// Scale every element of `v` by `s` in-place.
#[inline]
fn scale_vec_f32(v: &mut [f32], s: f32) {
    for x in v.iter_mut() {
        *x *= s;
    }
}

// ---------------------------------------------------------------------------
// Tests
// ---------------------------------------------------------------------------

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

    /// Helper: build a diagonally dominant tridiagonal matrix.
    fn tridiag_f32(n: usize, diag_val: f32, off_val: f32) -> CsrMatrix<f32> {
        let mut entries = Vec::new();
        for i in 0..n {
            entries.push((i, i, diag_val));
            if i > 0 {
                entries.push((i, i - 1, off_val));
            }
            if i + 1 < n {
                entries.push((i, i + 1, off_val));
            }
        }
        CsrMatrix::<f32>::from_coo(n, n, entries)
    }

    /// Helper: build a 3x3 system whose eigenvalues are in (0, 2) so that
    /// the Neumann series converges (rho(I - A) < 1).
    fn test_matrix_f64() -> CsrMatrix<f64> {
        CsrMatrix::<f64>::from_coo(
            3,
            3,
            vec![
                (0, 0, 1.0),
                (0, 1, -0.1),
                (1, 0, -0.1),
                (1, 1, 1.0),
                (1, 2, -0.1),
                (2, 1, -0.1),
                (2, 2, 1.0),
            ],
        )
    }

    #[test]
    fn test_new() {
        let solver = NeumannSolver::new(1e-8, 100);
        assert_eq!(solver.tolerance, 1e-8);
        assert_eq!(solver.max_iterations, 100);
    }

    #[test]
    fn test_spectral_radius_identity() {
        let identity = CsrMatrix::<f32>::identity(4);
        let rho = NeumannSolver::estimate_spectral_radius(&identity);
        assert!(rho < 0.1, "expected rho ~ 0 for identity, got {rho}");
    }

    #[test]
    fn test_spectral_radius_pure_diagonal() {
        // For a pure diagonal matrix D, D^{-1}A = I, so M = I - I = 0.
        // The spectral radius should be ~0.
        let a = CsrMatrix::<f32>::from_coo(
            3, 3,
            vec![(0, 0, 0.5_f32), (1, 1, 0.5), (2, 2, 0.5)],
        );
        let rho = NeumannSolver::estimate_spectral_radius(&a);
        assert!(rho < 0.1, "expected rho ~ 0 for diagonal matrix, got {rho}");
    }

    #[test]
    fn test_spectral_radius_empty() {
        let empty = CsrMatrix::<f32> {
            row_ptr: vec![0], col_indices: vec![], values: vec![],
            rows: 0, cols: 0,
        };
        assert_eq!(NeumannSolver::estimate_spectral_radius(&empty), 0.0);
    }

    #[test]
    fn test_spectral_radius_non_diag_dominant() {
        // Matrix where off-diagonal entries dominate:
        // [1  2]
        // [2  1]
        // D^{-1}A = [[1, 2], [2, 1]], so M = I - D^{-1}A = [[0, -2], [-2, 0]].
        // Eigenvalues of M are +2 and -2, so rho(M) = 2 > 1.
        let a = CsrMatrix::<f32>::from_coo(
            2, 2,
            vec![(0, 0, 1.0_f32), (0, 1, 2.0), (1, 0, 2.0), (1, 1, 1.0)],
        );
        let rho = NeumannSolver::estimate_spectral_radius(&a);
        assert!(rho > 1.0, "expected rho > 1 for non-diag-dominant matrix, got {rho}");
    }

    #[test]
    fn test_solve_identity() {
        let identity = CsrMatrix::<f32>::identity(3);
        let rhs = vec![1.0_f32, 2.0, 3.0];
        let solver = NeumannSolver::new(1e-6, 100);
        let result = solver.solve(&identity, &rhs).unwrap();
        for (i, (&e, &a)) in rhs.iter().zip(result.solution.iter()).enumerate() {
            assert!((e - a).abs() < 1e-4, "index {i}: expected {e}, got {a}");
        }
        assert!(result.residual_norm < 1e-6);
    }

    #[test]
    fn test_solve_diagonal() {
        let a = CsrMatrix::<f32>::from_coo(
            3, 3,
            vec![(0, 0, 0.5_f32), (1, 1, 0.5), (2, 2, 0.5)],
        );
        let rhs = vec![1.0_f32, 1.0, 1.0];
        let solver = NeumannSolver::new(1e-6, 200);
        let result = solver.solve(&a, &rhs).unwrap();
        for (i, &val) in result.solution.iter().enumerate() {
            assert!((val - 2.0).abs() < 0.01, "index {i}: expected ~2.0, got {val}");
        }
    }

    #[test]
    fn test_solve_tridiagonal() {
        // diag=1.0, off=-0.1: Jacobi iteration matrix has rho ~ 0.17.
        // Use 1e-6 tolerance since f32 accumulation limits floor.
        let a = tridiag_f32(5, 1.0, -0.1);
        let rhs = vec![1.0_f32, 0.0, 1.0, 0.0, 1.0];
        let solver = NeumannSolver::new(1e-6, 1000);
        let result = solver.solve(&a, &rhs).unwrap();
        assert!(result.residual_norm < 1e-4);
        assert!(result.iterations > 0);
        assert!(!result.convergence_history.is_empty());
    }

    #[test]
    fn test_solve_empty_system() {
        let a = CsrMatrix::<f32> {
            row_ptr: vec![0], col_indices: vec![], values: vec![],
            rows: 0, cols: 0,
        };
        let result = NeumannSolver::new(1e-6, 10).solve(&a, &[]).unwrap();
        assert_eq!(result.iterations, 0);
        assert!(result.solution.is_empty());
    }

    #[test]
    fn test_solve_dimension_mismatch() {
        let a = CsrMatrix::<f32>::identity(3);
        let rhs = vec![1.0_f32, 2.0];
        let err = NeumannSolver::new(1e-6, 100).solve(&a, &rhs).unwrap_err();
        let msg = err.to_string();
        assert!(msg.contains("dimension") || msg.contains("mismatch"), "got: {msg}");
    }

    #[test]
    fn test_solve_non_square() {
        let a = CsrMatrix::<f32>::from_coo(2, 3, vec![(0, 0, 1.0_f32), (1, 1, 1.0)]);
        let rhs = vec![1.0_f32, 1.0];
        let err = NeumannSolver::new(1e-6, 100).solve(&a, &rhs).unwrap_err();
        let msg = err.to_string();
        assert!(msg.contains("square") || msg.contains("dimension"), "got: {msg}");
    }

    #[test]
    fn test_solve_divergent_matrix() {
        // Non-diag-dominant: off-diagonal entries larger than diagonal.
        let a = CsrMatrix::<f32>::from_coo(
            2, 2,
            vec![(0, 0, 1.0_f32), (0, 1, 2.0), (1, 0, 2.0), (1, 1, 1.0)],
        );
        let rhs = vec![1.0_f32, 1.0];
        let err = NeumannSolver::new(1e-6, 100).solve(&a, &rhs).unwrap_err();
        assert!(err.to_string().contains("spectral radius"), "got: {}", err);
    }

    #[test]
    fn test_convergence_history_monotone() {
        let a = CsrMatrix::<f32>::identity(4);
        let rhs = vec![1.0_f32; 4];
        let result = NeumannSolver::new(1e-10, 50).solve(&a, &rhs).unwrap();
        assert!(!result.convergence_history.is_empty());
        for window in result.convergence_history.windows(2) {
            assert!(
                window[1].residual_norm <= window[0].residual_norm + 1e-12,
                "residual not decreasing: {} -> {}",
                window[0].residual_norm, window[1].residual_norm,
            );
        }
    }

    #[test]
    fn test_algorithm_tag() {
        let a = CsrMatrix::<f32>::identity(2);
        let rhs = vec![1.0_f32; 2];
        let result = NeumannSolver::new(1e-6, 100).solve(&a, &rhs).unwrap();
        assert_eq!(result.algorithm, Algorithm::Neumann);
    }

    #[test]
    fn test_solver_engine_trait_f64() {
        let solver = NeumannSolver::new(1e-6, 200);
        let engine: &dyn SolverEngine = &solver;
        let a = test_matrix_f64();
        let rhs = vec![1.0_f64, 0.0, 1.0];
        let budget = ComputeBudget::default();
        let result = engine.solve(&a, &rhs, &budget).unwrap();
        assert!(result.residual_norm < 1e-4);
        assert_eq!(result.algorithm, Algorithm::Neumann);
    }

    #[test]
    fn test_larger_system_accuracy() {
        let n = 50;
        // diag=1.0, off=-0.1: Jacobi-preconditioned Neumann converges.
        // Use 1e-6 tolerance for f32 precision headroom.
        let a = tridiag_f32(n, 1.0, -0.1);
        let rhs: Vec<f32> = (0..n).map(|i| (i as f32 + 1.0) / n as f32).collect();
        let result = NeumannSolver::new(1e-6, 2000).solve(&a, &rhs).unwrap();
        assert!(result.residual_norm < 1e-6, "residual too large: {}", result.residual_norm);
        let mut ax = vec![0.0f32; n];
        a.spmv(&result.solution, &mut ax);
        for i in 0..n {
            assert!((ax[i] - rhs[i]).abs() < 1e-4, "A*x[{i}]={} but b[{i}]={}", ax[i], rhs[i]);
        }
    }

    #[test]
    fn test_scalar_system() {
        let a = CsrMatrix::<f32>::from_coo(1, 1, vec![(0, 0, 0.5_f32)]);
        let rhs = vec![4.0_f32];
        let result = NeumannSolver::new(1e-8, 200).solve(&a, &rhs).unwrap();
        assert!((result.solution[0] - 8.0).abs() < 0.01, "expected ~8.0, got {}", result.solution[0]);
    }

    #[test]
    fn test_estimate_complexity() {
        let solver = NeumannSolver::new(1e-6, 1000);
        let profile = SparsityProfile {
            rows: 100, cols: 100, nnz: 500, density: 0.05,
            is_diag_dominant: true, estimated_spectral_radius: 0.5,
            estimated_condition: 3.0, is_symmetric_structure: true,
            avg_nnz_per_row: 5.0, max_nnz_per_row: 8,
        };
        let estimate = solver.estimate_complexity(&profile, 100);
        assert_eq!(estimate.algorithm, Algorithm::Neumann);
        assert!(estimate.estimated_flops > 0);
        assert!(estimate.estimated_iterations > 0);
        assert!(estimate.estimated_memory_bytes > 0);
        assert_eq!(estimate.complexity_class, ComplexityClass::SublinearNnz);
    }
}