Module neumann 

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:

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:

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.

Structs

NeumannSolver

Neumann Series solver for sparse linear systems.