Repository navigation
Several bugs in iterative solvers #603
Description
Several of the iterative solvers in psydac.linalg.solvers crash on edge cases, or return wrong results for complex operators. In one case the solver reports convergence for an unconverged solution. All of the bugs below were reproduced with NumPy 2.5.
| Solver | Bug | Symptom |
|---|---|---|
MinimumResidual |
inf is used but never imported |
NameError when the operator norm or the solution norm vanishes (e.g. zero operator) |
LSMR |
np.infty was removed in NumPy 2.0 (we require numpy >= 2.1) |
AttributeError when the residual vanishes exactly (e.g. A = 2*I, b = 1) |
GMRES |
Off-by-one: the solution is built from the first k Arnoldi vectors only, dropping the last one |
The returned x is always one iteration behind the reported residual. With maxiter = n, success=True and res_norm ≈ 1e-15 are reported, while the true residual is ≈ 1e-4 |
GMRES |
Wrong for complex operators: the Gram-Schmidt coefficients are p.inner(Q[i]) instead of Q[i].inner(p), and the Givens rotations are real-only |
The Krylov basis is not orthogonal. A complex 12×12 system needs 25 iterations instead of at most 12 |
BiConjugateGradientStabilized (with preconditioner) |
The inner products with the shadow residual rp0 conjugate the wrong vector (vp.inner(rp0), rp.inner(rp0)) |
Fails for complex operators (residual ≈ 7e-2 after 100 iterations). This is why test_solver_tridiagonal skipped this case |
ConjugateGradient (with and without preconditioner) |
The loop counter starts at 2 and is never assigned when maxiter=1 |
UnboundLocalError |
LSMR |
The port from SciPy dropped the early exits for b = 0 and A^H r0 = 0 |
ZeroDivisionError |
Why the tests did not catch these
test_solver_tridiagonaluses a generousmaxiterand only checks the final error, so the GMRES off-by-one never shows up.- Complex operators are skipped for preconditioned BiCGSTAB, and the complex GMRES test passes only because
maxiteris large. - No test covers
maxiter=1,b = 0or degenerate operators.
Out of scope
These are related but not addressed here:
MinimumResidualreportssuccess=Truefor an unsolvable system (zero operator), becausesuccessrelies on the internal residual estimate, which vanishes at a Lanczos breakdown.- BiCG and LSMR fail for complex operators (see LSMR and BiCG currently failing tests for complex matrices #557).
- The meaning of
niterdiffers between solvers: CG and GMRES count the initial residual as iteration 1, while the others count only the iterations performed (see Inconsistent meaning ofniterandmaxiterin iterative solvers #630)
Activity
Remove GMRES success assertion, which doesn't test off-by-one bugfix (#…