Fix several bugs in iterative linear solvers by yguclu · Pull Request #604 · pyccel/psydac

yguclu · 2026-10-05T18:18:42Z

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

MINRES and LSMR: import inf from math, which MINRES used without importing. Replace np.infty, removed in NumPy 2.0, in LSMR.
GMRES: build the solution from all k+1 Arnoldi vectors instead of k. The returned x now matches the reported residual.
GMRES: support complex operators by using the Gram-Schmidt coefficients Q[i]^H p and unitary Givens rotations.
Preconditioned BiCGSTAB: conjugate the shadow residual rp0 in the inner products, as in the unpreconditioned version. This fixes complex operators.
CG: initialize the loop counter before the loop, so that maxiter=1 no longer raises UnboundLocalError.
LSMR: restore SciPy's early exits for b = 0 (return x = 0) and A^H (b - A x0) = 0 (return x0), which previously caused a ZeroDivisionError.

Behavior notes

niter is unchanged for CG and GMRES. Both still count the initial residual as iteration 1. Because of this, GMRES reports niter = maxiter + 1 when it does not converge.
LSMR can now return without iterating (niter = 0) when b = 0 or when x0 is 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 with 2*I
test_GMRES_solve: compares the reported residual with the true residual, for real and complex operators
test_ConjugateGradient_solve_maxiter_1: CG with and without preconditioner
test_LSMR_solve_early_exit

The 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:

MINRES reports success=True for an unsolvable system.
BiCG and LSMR still fail for complex operators (LSMR and BiCG currently failing tests for complex matrices #557).
The meaning of niter differs between solvers.

codacy-production · 2026-10-05T18:19:31Z

Up to standards ✅

🟢 Issues 0 issues

Results:
0 new issues

View in Codacy

🟢 Metrics 4 complexity · 0 duplication

Metric Results

Complexity 4

Duplication 0

View in Codacy

🟢 Coverage 97.37% diff coverage · +0.01% coverage variation

Metric Results

Coverage variation ✅ +0.01% coverage variation

Diff coverage ✅ 97.37% diff coverage

View coverage diff in Codacy

Coverage variation details

Coverable lines Covered lines Coverage

Common ancestor commit (f2da0cd) 51844 33084 63.81%

Head commit (778f084) 103790 (+51946) 66242 (+33158) 63.82% (+0.01%)

Coverage variation is the difference between the coverage for the head and common ancestor commits of the pull request branch: <coverage of head commit> - <coverage of common ancestor commit>

Diff coverage details

Coverable lines Covered lines Diff coverage

Pull request (#604) 76 74 97.37%

Diff coverage is the percentage of lines that are covered by tests out of the coverable lines that the pull request added or modified: <covered lines added or modified>/<coverable lines added or modified> * 100%

_{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.}

`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`.

campospinto

thanks for this PR which seems good for me.

I have mostly one question:

the PR descriptions states that it fixes #603 but also says that a few issues listed in #603 are not addressed in this PR (eg MINRES reports success=True for an unsolvable system), so there should probably be two different issues, one fixed and one not yet?

jowezarek

I added a few comments related to the new tests (and a remark related to define_data). My main concern is that LSMR returning success in case of an infeasible system (0 matrix) is misleading. I did not check the solver details as I am not familiar with the theory.

jowezarek · 2026-10-16T08:17:23Z

    s = V.starts[0]

    # Build banded matrix with 2p+1 diagonals: must be symmetric and positive definite
    # Here we assign value 2*p on main diagonal and -1 on other diagonals


I suppose this function would benefit from updated comments / docstrings.
It's no longer true that the main diagonal entries are 2p and the off-diagonal entries are -1. Now, one can pass matrix_data (length >= 3), but only the first 3 values are used to fill the matrix. Also, the resulting matrix could be not symmetric and not positive definit.

jowezarek · 2026-10-16T08:34:39Z

+
+    r = b - A @ x
+    assert np.isclose(info['res_norm'], np.sqrt(r.inner(r).real), rtol=1e-8, atol=1e-12)
+    assert info['success'] == (maxiter == n)


I understand the previous assert (test whether the true residual is reported in info['res_norm'], but what is the purpose of this assertion? A and b are chosen in such a way that the solver requires >= 12 iterations (1+12=13 niter) to converge? But why is this checked?

jowezarek · 2026-10-16T08:42:24Z

+
+#===============================================================================
+# With maxiter=1, CG only evaluates the residual of the initial guess
+@pytest.mark.parametrize('use_jacobi_pc', [False, True])


Regarding Behavior notes of the PR description: GMRES also reports niter=maxiter+1 in case of convergence in case of convergence with the last iteration (see test above). Question: Why not change it such that niter=1 actually does exactly 1 iteration (from initial guess to improved guess)?

jowezarek · 2026-10-16T08:54:11Z

+    expected = xe if rhs else V.zeros()
+    assert np.array_equal(x.toarray(), expected.toarray())
+    assert info['niter'] == 0
+    assert info['success']


Now LSMR reports success in case of a zero matrix and a non-zero rhs. Is that really intended?

e-moral-sanchez

I compared the changes against PETSc implementations of BiCGStab and GMRES and found some discrepancies.

e-moral-sanchez · 2026-10-16T09:05:20Z

+
+    x = inverse(A, solver, tol=1e-10) @ b
+
+    assert np.array_equal(x.toarray(), np.full(6, expected))


Why would the expected solution value depend on the solver if it converges? I don't think this test brings anything.

e-moral-sanchez · 2026-10-16T09:09:34Z

+
+    r = b - A @ x
+    assert np.isclose(info['res_norm'], np.sqrt(r.inner(r).real), rtol=1e-8, atol=1e-12)
+    assert info['success'] == (maxiter == n)


This is only true with exact precision. Also here we set tol=1e-10. I don't think this test brings anything.

e-moral-sanchez · 2026-10-16T09:10:55Z

+
+    assert np.array_equal(x.toarray(), V.zeros().toarray())
+    assert info['niter'] == 1
+    assert not info['success']


Do we need to test this?

e-moral-sanchez · 2026-10-16T09:11:42Z

+    assert np.array_equal(x.toarray(), expected.toarray())
+    assert info['niter'] == 0
+    assert info['success']
+


I also would not add this test, it brings nothing.

e-moral-sanchez · 2026-10-16T09:42:43Z

            A.dot(pp, out=v)
            pc.dot(v, out=vp)
-            alphap = rhop / vp.inner(rp0)
+            alphap = rhop / rp0.inner(vp)


I agree with this change, but I think it's worth commenting it. In PSYDAC inner takes the conjugate of the first argument (

psydac/psydac/linalg/kernels/inner_kernels.py

Line 39 in daec6f7

res += v1[i0].conjugate() * v2[i0]

) . In PETSc, they take it in the second argument (https://petsc.org/release/manualpages/Vec/VecDot/) and in their implementation they do "VecDot(V, RP, &d1))" see line 78 in https://petsc.org/release/src/ksp/ksp/impls/bcgs/bcgs.c.html#KSPBCGS

e-moral-sanchez · 2026-10-16T11:21:31Z

            test1 = normr / normb
            if (normA * normr) != 0:test2 = normar / (normA * normr)
-            else:test2 = np.infty
+            else:test2 = inf


Is this change necessary?

e-moral-sanchez · 2026-10-16T11:29:30Z

            # update the residual vector
            beta.append(- sn[k] * beta[k])
-            beta[k] *= cn[k]
+            beta[k] *= cn[k].conjugate()


I agree with this change, and here also I think it's good to comment that it is consistent with the PETSc implementation: see line 377 of https://petsc.org/release/src/ksp/ksp/impls/gmres/gmres.c.html#KSPGMRES

e-moral-sanchez · 2026-10-16T11:42:09Z

-        h[k] *= cn[k]
-        h[k] += sn[k] * h[k+1]
+        h[k] *= cn[k].conjugate()
+        h[k] += sn[k].conjugate() * h[k+1]


This one is not conjugated in PETSc GMRES: see line 378 in https://petsc.org/release/src/ksp/ksp/impls/gmres/gmres.c.html#KSPGMRES

e-moral-sanchez · 2026-10-16T11:45:25Z

-            h[i] *= cn[i]
-            h[i] += sn[i] * h[i+1]
+            h[i] *= cn[i].conjugate()
+            h[i] += sn[i].conjugate() * h[i+1]


This one is not conjugated in PETSc GMRES: see line 378 in https://petsc.org/release/src/ksp/ksp/impls/gmres/gmres.c.html#KSPGMRES

yguclu changed the title ~~Fix several bugs in linear solvers~~ Fix several bugs in iterative linear solvers Oct 5, 2026

yguclu linked an issue Oct 6, 2026 that may be closed by this pull request

Several bugs in iterative solvers #603

Open

yguclu force-pushed the yguclu-fix-linear-solvers branch from 7cae052 to caf0703 Compare October 6, 2026 08:31

yguclu added 6 commits October 9, 2026 21:52

yguclu force-pushed the yguclu-fix-linear-solvers branch from caf0703 to dd703ea Compare October 9, 2026 20:00

yguclu changed the base branch from yguclu-improve-Geometry to devel October 9, 2026 20:01

yguclu requested review from a team, e-moral-sanchez and jowezarek October 9, 2026 21:58

yguclu marked this pull request as ready for review October 9, 2026 21:59

yguclu added the Next Release Must be in next release label Oct 9, 2026

Mention fixes of iterative linear solvers in CHANGELOG.md

778f084

yguclu marked this pull request as draft October 9, 2026 22:23

yguclu marked this pull request as ready for review October 9, 2026 22:23

yguclu removed the request for review from a team October 15, 2026 12:36

yguclu self-assigned this Oct 15, 2026

campospinto reviewed Oct 15, 2026

View reviewed changes

Comment thread psydac/linalg/solvers.py

Comment thread psydac/linalg/solvers.py

jowezarek reviewed Oct 16, 2026

View reviewed changes

e-moral-sanchez requested changes Oct 16, 2026

View reviewed changes


		x = inverse(A, solver, tol=1e-10) @ b

		assert np.array_equal(x.toarray(), np.full(6, expected))

Conversation

yguclu commented Oct 5, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Changes

Behavior notes

Tests

Not addressed

Uh oh!

codacy-production Bot commented Oct 5, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Up to standards ✅

Uh oh!

campospinto left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

jowezarek left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

e-moral-sanchez left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

This branch has not been deployed

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

yguclu commented Oct 5, 2026 •

edited

Loading

codacy-production Bot commented Oct 5, 2026 •

edited

Loading