Skip to main content

Module poly

Module poly 

Source
Expand description

Polynomial preconditioner (Neumann series and Chebyshev).

A polynomial preconditioner approximates A^{-1} by a polynomial in A: M^{-1} = p(A). Applying it is nothing but a handful of sparse matrix-vector products and vector updates — there are no triangular solves. That makes it the odd one out in this crate: where ILU/IC apply a sequential forward/back substitution, a polynomial preconditioner is built entirely from A * x, which parallelises and vectorises freely.

Two flavours are provided:

  • Neumann seriesM^{-1} = w * sum_{k=0}^{d} (I - w A)^k. One real damping parameter w; converges when 0 < w * lambda < 2 across the spectrum.
  • Chebyshev — the degree-d Chebyshev polynomial that minimises max |1 - lambda p(lambda)| over [lambda_min, lambda_max]. Sharper than Neumann for the same degree, but it needs an estimate of the spectral interval.

§When to use it

Reach for a polynomial preconditioner when matrix-vector products are cheap and plentiful but triangular solves are a bottleneck — many cores, a GPU, or a distributed operator where the sequential sweep of an ILU does not scale. It is also the classic choice for a smoother inside multigrid. On a single core it rarely beats crate::Ic0 / crate::Ilu0; its appeal is parallelism and the absence of any factorisation.

Chebyshev assumes a Hermitian positive-definite operator and is only as good as its [lambda_min, lambda_max] estimate: an over-estimated lambda_min (or under-estimated lambda_max) degrades or even diverges the polynomial. Pass BoundEstimate::Manual when you know the spectrum; otherwise BoundEstimate::PowerIteration gives a tight lambda_max and a conservative lambda_min.

§Storage

The operator A is stored as an owned CSC copy (apply reads it through a faer::sparse::SparseColMatRef). The recurrence’s temporaries — one work column for Neumann, two for Chebyshev — come from the caller’s MemStack, so apply allocates no heap memory.

Structs§

Poly
Polynomial preconditioner M^{-1} = p(A).
PolyParams
Tuning parameters for Poly::try_new.

Enums§

BoundEstimate
How to obtain the spectral interval for Poly::try_new_auto.
PolyError
Error returned by polynomial-preconditioner construction.
PolyKind
Which polynomial to use for M^{-1} = p(A).