Repository navigation
Conversation
Up to standards ✅🟢 Issues
|
| Metric | Results |
|---|---|
| Complexity | 4 |
| Duplication | 0 |
NEW Get contextual insights on your PRs based on Codacy's metrics, along with PR and Jira context, without leaving GitHub. Enable AI reviewer
TIP This summary will be updated as you push new changes.
yguclu
changed the title
`MinimumResidual.solve` used `inf` without importing it, raising a NameError when the operator norm or the solution norm vanishes. `LSMR.solve` used `np.infty`, which was removed in NumPy 2.0, raising an AttributeError when the residual vanishes exactly. Both now use `math.inf`. Add `test_solver_diagonal`, which reaches these branches with diagonal operators (zero operator for MINRES, 2*I for LSMR).
`GMRES.solve` built the solution from the first k Arnoldi vectors only, dropping the last one. The returned x was therefore always one iteration behind the reported residual, and `success` could be True for an unconverged solution. Check convergence at the end of each iteration and use all k+1 Arnoldi vectors. The `niter` convention is unchanged. Add `test_GMRES_solve`, which compares the reported residual with the true residual b - A x.
The Arnoldi process computed the Gram-Schmidt coefficients as `p.inner(Q[i])`, i.e. conj(p) . Q[i], instead of `Q[i].inner(p)`, so the Krylov basis was not orthogonal for complex operators. The Givens rotations were also written for real numbers only. Use the coefficients Q[i]^H p and the unitary rotation [[conj(c), conj(s)], [-s, c]]. Complex systems now converge in at most n iterations, as in the real case. Extend `test_GMRES_solve` to complex operators.
`BiConjugateGradientStabilized.solve_with_pc` computed the inner products with the shadow residual rp0 as `vp.inner(rp0)` and `rp.inner(rp0)`, conjugating the wrong vector. This conjugated alpha and beta, so the solver failed for complex operators. Use `rp0.inner(vp)` and `rp0.inner(rp)`, as in the unpreconditioned version. Remove the corresponding skip from `test_solver_tridiagonal`.
In `ConjugateGradient.solve_without_pc` and `solve_with_pc` the loop counter starts at 2, because iteration 1 is the initial residual. With maxiter=1 the loop body never runs, and the counter was used unassigned to fill `niter`, raising an UnboundLocalError. Initialize the counter to 1 before the loop; the `niter` convention is unchanged. Add `test_ConjugateGradient_solve_maxiter_1`.
The SciPy port of `LSMR.solve` dropped the early exits taken before the main loop. With b = 0 the stopping test divided by norm(b) = 0, and when A^H (b - A x0) = 0 the plane rotations divided by zero. Now return x = 0 if b = 0, and return x0 if it is already a least-squares solution. Add `test_LSMR_solve_early_exit`.
yguclu
force-pushed
the
yguclu-fix-linear-solvers
branch
from
October 6, 2026 08:31
7cae052 to
caf0703
Compare
This branch has not been deployed
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Fix several bugs in the iterative solvers of
psydac.linalg.solvers, as reported in #603. Each bug is fixed in a separate commit, together with a regression test that fails before the fix.Fixes #603
Changes
inffrommath, which MINRES used without importing. Replacenp.infty, removed in NumPy 2.0, in LSMR.k+1Arnoldi vectors instead ofk. The returnedxnow matches the reported residual.Q[i]^H pand unitary Givens rotations.rp0in the inner products, as in the unpreconditioned version. This fixes complex operators.maxiter=1no longer raisesUnboundLocalError.b = 0(returnx = 0) andA^H (b - A x0) = 0(returnx0), which previously caused aZeroDivisionError.Behavior notes
niteris unchanged for CG and GMRES. Both still count the initial residual as iteration 1. Because of this, GMRES reportsniter = maxiter + 1when it does not converge.niter = 0) whenb = 0or whenx0is already a least-squares solution.Tests
New tests in
psydac/linalg/tests/test_solvers.py:test_solver_diagonal: MINRES with a zero operator, and LSMR with2*Itest_GMRES_solve: compares the reported residual with the true residual, for real and complex operatorstest_ConjugateGradient_solve_maxiter_1: CG with and without preconditionertest_LSMR_solve_early_exitThe skip of complex operators for preconditioned BiCGSTAB was removed from
test_solver_tridiagonal.Not addressed
These are listed in #603 and left for follow-up work:
success=Truefor an unsolvable system.niterdiffers between solvers.Notes
PR can be merged only after #527!