From 620b1c478c4deea232197c7d2cad9dc5c87fa32a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Mon, 5 Oct 2026 12:47:28 +0200 Subject: [PATCH 01/11] Add CLAUDE.md configuration file --- CLAUDE.md | 74 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 74 insertions(+) create mode 100644 CLAUDE.md diff --git a/CLAUDE.md b/CLAUDE.md new file mode 100644 index 000000000..d25c1bb03 --- /dev/null +++ b/CLAUDE.md @@ -0,0 +1,74 @@ +# Project: PSYDAC + +## Overview +Python 3 library for isogeometric analysis (IGA). Solve general systems of partial +differential equations (PDEs) in weak form, defined using the domain-specific +language provided by SymPDE. Supports finite element exterior calculus (FEEC) +with tensor-product spline spaces. Handles multi-patch geometries in various ways; +usually broken-FEEC a.k.a. CONGA (conforming/non-conforming Galerkin). +Python code is automatically generated for the assembly of user-defined functionals, +linear forms, and bilinear forms. This Python code is then accelerated to C/Fortran +speed using Pyccel. The library enables large parallel computations on distributed- +memory supercomputers using MPI and OpenMP. + +## Tech stack +- Python 3.10+, type hints required everywhere +- meson-python + mesonpy for building +- pip with submodules (for igakit) +- pytest + pytest-cov + pytest-mpi + pytest-xdist for testing +- sympde for symbolic definition of weak formulations +- pyccel for translation of Python kernels to Fortran or C +- mpi4py for MPI parallelization +- h5py for parallel I/O +- petsc4py for direct linear solvers (optional install) + +## Code standards +- Follow PEP 8 style guide whenever possible +- Docstrings of public functions and classes follow Numpydoc conventions +- pylint + black + isort for linting/typing +- Computational kernels to be accelerated with Pyccel named as `_kernels.py` +- Minimize code duplication. Never add code that already exists +- Strive for concise, clean, and human-readable code. Avoid useless verbosity +- Self-explanatory names for variables, functions, and classes +- Do not reassign value to an existing variable, especially if the type changes +- Use short comments (one-liners or inline) to explain "why" rather than "what" or "how" + +## Testing conventions +- Each library subpackage contains a `tests/` folder with an `__init__.py` file +- Unit tests are part of the library and shipped with it +- Unit tests can be run with `psydac test` CLI command (see `README.md`) +- Test file names follow `test_.py` +- Test function names follow `test_` or `test__` +- Parametrize unit tests with `@pytest.parametrize` to minimize code duplication +- Keep run time of tests at a minimum +- Aim for 100% coverage on newly committed code + +## File structure +- psydac/api - high-level Python interface +- psydac/cad - computer-aided design (CAD) functionality +- psydac/cmd - command-line interface (CLI) commands +- psydac/core - splines functionality +- psydac/ddm - MPI decomposition of domain, vectors of spline coefficients, and matrices +- psydac/feec - finite element exterior calculus (FEEC) +- psydac/fem - "finite element method" middle-level Python interface (spaces and fields) +- psydac/linalg - "linear algebra" low-level Python interface (vectors, linear operators, iterative solvers) +- psydac/polar - implementation of polar splines for H1 spaces +- psydac/pyccel - distillation of an old version of Pyccel for Python code generation (to be removed) +- psydac/utilities - various generic functions used across the library +- docs - Sphinx documentation +- examples - Various examples of library usage (mostly Jupyter notebooks) +- scripts - Various scripts for library maintenance +- subprojects - Git submodules + +## Always do +- Before claiming done: + . Run `python -m black --check $(git ls-files "*.py")` + . Run `python -m isort --check $(git ls-files "*.py" | grep -v "FOUND_DUPLICATED_IMPORT.py")` + . Run `pylint --disable=all --enable=unused-import ./psydac --ignore-paths="__pyccel*" --ignore-paths="__epyccel*"` +- Unless specified differently, assume process with `rank=0` in MPI communicator is root +- Use root MPI process to perform serial operations (e.g. open/save image) + +## Never do +- Commit secrets or `.env` files +- Add new dependencies without team discussion +- Use `eval` function together with `input` From e37983cb1e95ce57874d0245d571be7a6d8ec3d9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Mon, 5 Oct 2026 17:33:38 +0200 Subject: [PATCH 02/11] Remove blank imports from psydac/cad/__init__.py --- psydac/cad/__init__.py | 7 ------- 1 file changed, 7 deletions(-) diff --git a/psydac/cad/__init__.py b/psydac/cad/__init__.py index 6f39d9a03..419109b64 100644 --- a/psydac/cad/__init__.py +++ b/psydac/cad/__init__.py @@ -3,10 +3,3 @@ # LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE # # for full license details. # #---------------------------------------------------------------------------# - -__all__ = ['geometry'] - -from psydac.cad import geometry -from psydac.cad import cad -from psydac.cad import gallery -from psydac.cad import utils From 03f1d32e4c1377bbdf9de3ef73ae86f9d2a3c43c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Mon, 5 Oct 2026 17:37:18 +0200 Subject: [PATCH 03/11] Remove obsolete test_geometry_1() from psydac.cad.tests.test_geometry --- psydac/cad/tests/test_geometry.py | 8 -------- 1 file changed, 8 deletions(-) diff --git a/psydac/cad/tests/test_geometry.py b/psydac/cad/tests/test_geometry.py index 18ad343b1..1811fde14 100644 --- a/psydac/cad/tests/test_geometry.py +++ b/psydac/cad/tests/test_geometry.py @@ -365,14 +365,6 @@ def test_import_geopdes_to_nurbs(ncells, degree): if isinstance(mapping, NurbsMapping): assert np.allclose(L_shaped.weights.flatten(), mapping._weights_field.coeffs.toarray(), 1e-15, 1e-15) -#============================================================================== -@pytest.mark.xfail -def test_geometry_1(): - - line = Geometry.as_line(ncells=[10]) - square = Geometry.as_square(ncells=[10, 10]) - cube = Geometry.as_cube(ncells=[10, 10, 10]) - #============================================================================== # CLEAN UP SYMPY NAMESPACE #============================================================================== From 68bb98efea95f78a598cc3f300d9a868b61de0e1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Mon, 5 Oct 2026 18:24:28 +0200 Subject: [PATCH 04/11] Add Pylint configuration to pyproject.toml --- pyproject.toml | 34 ++++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) diff --git a/pyproject.toml b/pyproject.toml index ad867dc80..36c290a3e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -76,3 +76,37 @@ psydac = "psydac.cmd.main:psydac_command" setup = ['--default-library=static'] dist = ['--include-subprojects'] install = ['--only=python-modules'] + +[tool.pylint.main] +py-version = '3.10' +ignore-paths = ['psydac/pyccel/.*', '.*/__e?pyccel__/.*', '.*/__psydac__/.*'] +ignored-modules = ['mpi4py'] # C extension, members not visible to Pylint + +[tool.pylint.basic] +no-docstring-rgx = '^(_|test_|Test)' + +[tool.pylint.'messages control'] +disable = [ + # formatting owned by black / isort + 'line-too-long', + 'trailing-whitespace', + 'bad-indentation', + 'multiple-statements', + 'missing-final-newline', + 'trailing-newlines', + 'wrong-import-order', + # style choices for mathematical code + 'invalid-name', + 'too-many-arguments', + 'too-many-positional-arguments', + 'too-many-locals', + 'too-many-branches', + 'too-many-statements', + 'too-many-instance-attributes', + 'too-few-public-methods', + 'too-many-nested-blocks', + # pytest fixtures shadow outer names by design + 'redefined-outer-name', + 'duplicate-code', + 'fixme', +] From 8730761234a843ef389e58aea3b795cfdf9d54a3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Mon, 5 Oct 2026 18:28:01 +0200 Subject: [PATCH 05/11] Update CHANGELOG.md --- CHANGELOG.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index ad4d4aba0..f9dc0e1c3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -14,6 +14,8 @@ All notable changes to this project will be documented in this file. - #567 : Improve `psydac test` command (with several new features) - #565 : Expand editable install info in `README.md` - [DEVELOPER] Create action `install_petsc4py` to install PETSc & `petsc4py` w/ complex support +- [DEVELOPER] Configure Pylint in `pyproject.toml` +- [DEVELOPER] Add `CLAUDE.md` with developer rules ### Fixed From bf11e434eb6952fd1118ccdfbfcf4db77d375e56 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Mon, 5 Oct 2026 19:06:21 +0200 Subject: [PATCH 06/11] Fix undefined `inf` in MINRES and remove `np.infty` in LSMR `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). --- psydac/linalg/solvers.py | 19 ++++++++++++------- psydac/linalg/tests/test_solvers.py | 15 +++++++++++++++ 2 files changed, 27 insertions(+), 7 deletions(-) diff --git a/psydac/linalg/solvers.py b/psydac/linalg/solvers.py index fc3cce3cc..b7d8c87fa 100644 --- a/psydac/linalg/solvers.py +++ b/psydac/linalg/solvers.py @@ -7,7 +7,7 @@ This module provides iterative solvers and preconditioners. """ -from math import sqrt +from math import inf, sqrt import numpy as np import warnings @@ -1218,12 +1218,17 @@ def solve(self, b, out=None): ynorm = sqrt(x.inner(x)) rnorm = phibar - if ynorm == 0 or Anorm == 0:test1 = inf - #else:test1 = rnorm / (Anorm*ynorm) # ||r|| / (||A|| ||x||) - else:test1 = rnorm # ||r|| + if ynorm == 0 or Anorm == 0: + test1 = inf + #else: + # test1 = rnorm / (Anorm*ynorm) # ||r|| / (||A|| ||x||) + else: + test1 = rnorm # ||r|| - if Anorm == 0:test2 = inf - else:test2 = root / Anorm # ||Ar|| / (||A|| ||r||) + if Anorm == 0: + test2 = inf + else: + test2 = root / Anorm # ||Ar|| / (||A|| ||r||) # Estimate cond(A). # In this version we look at the diagonals of R in the @@ -1596,7 +1601,7 @@ def solve(self, b, out=None): test1 = normr / normb if (normA * normr) != 0:test2 = normar / (normA * normr) - else:test2 = np.infty + else:test2 = inf test3 = 1 / condA t1 = test1 / (1 + normA * normx / normb) rtol = btol + atol * normA * normx / normb diff --git a/psydac/linalg/tests/test_solvers.py b/psydac/linalg/tests/test_solvers.py index ee122d876..2e0e50f16 100644 --- a/psydac/linalg/tests/test_solvers.py +++ b/psydac/linalg/tests/test_solvers.py @@ -252,6 +252,21 @@ def test_solver_tridiagonal(n, p, dtype, solver, use_jacobi_pc, verbose=False): assert errh_norm < tol assert (solver == 'CG' and use_jacobi_pc) or errc_norm < tol +#=============================================================================== +# Diagonal operators make the stopping tests divide by a vanishing norm +@pytest.mark.parametrize( + ('solver', 'diagonal', 'expected'), [('MINRES', 0.0, 0.0), ('LSMR', 2.0, 0.5)] +) +def test_solver_diagonal(solver: str, diagonal: float, expected: float) -> None: + + V, A, _ = define_data(6, 1, [0.0, diagonal, 0.0]) + b = V.zeros() + b[:] = 1.0 + + x = inverse(A, solver, tol=1e-10) @ b + + assert np.array_equal(x.toarray(), np.full(6, expected)) + #=============================================================================== def test_LST_preconditioner(comm=None): From 1d38353eece651e82c92c9bf6467228a7a83dc5a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Mon, 5 Oct 2026 19:20:24 +0200 Subject: [PATCH 07/11] Fix off-by-one in GMRES solution update `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. --- psydac/linalg/solvers.py | 15 ++++++++------- psydac/linalg/tests/test_solvers.py | 17 +++++++++++++++++ 2 files changed, 25 insertions(+), 7 deletions(-) diff --git a/psydac/linalg/solvers.py b/psydac/linalg/solvers.py index b7d8c87fa..82a459a5b 100644 --- a/psydac/linalg/solvers.py +++ b/psydac/linalg/solvers.py @@ -1782,8 +1782,6 @@ def solve(self, b, out=None): # Iterate to convergence for k in range(maxiter): - if am < tol: - break # run Arnoldi self.arnoldi(k, p) @@ -1799,16 +1797,19 @@ def solve(self, b, out=None): if verbose: print( template.format( k+2, am ) ) + if am < tol: + break + if verbose: - print( "+---------+---------------------+") - # calculate result - y = self.solve_triangular(self._H[:k, :k], beta[:k]) # system of upper triangular matrix + print( "+---------+---------------------+") + # calculate result from all k+1 Arnoldi vectors + y = self.solve_triangular(self._H[:k+1, :k+1], beta[:k+1]) # system of upper triangular matrix - for i in range(k): + for i in range(k+1): x.mul_iadd(y[i], self._Q[i]) # Convergence information - self._info = {'niter': k+1, 'success': bool(am < tol), 'res_norm': am} + self._info = {'niter': k+2, 'success': bool(am < tol), 'res_norm': am} if recycle: x.copy(out=self._options["x0"]) diff --git a/psydac/linalg/tests/test_solvers.py b/psydac/linalg/tests/test_solvers.py index 2e0e50f16..8edf38abc 100644 --- a/psydac/linalg/tests/test_solvers.py +++ b/psydac/linalg/tests/test_solvers.py @@ -267,6 +267,23 @@ def test_solver_diagonal(solver: str, diagonal: float, expected: float) -> None: assert np.array_equal(x.toarray(), np.full(6, expected)) +#=============================================================================== +# GMRES converges in at most n iterations, and must report the true residual +@pytest.mark.parametrize('maxiter', [1, 6, 12]) +def test_GMRES_solve(maxiter: int) -> None: + + n = 12 + _, A, xe = define_data(n, 1, [-7, -6, -1]) + b = A @ xe + + solver = inverse(A, 'GMRES', tol=1e-10, maxiter=maxiter) + x = solver @ b + info = solver.get_info() + + 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) + #=============================================================================== def test_LST_preconditioner(comm=None): From aba271c9ce93c2d600e42af35134147b4c2f4233 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Mon, 5 Oct 2026 19:26:19 +0200 Subject: [PATCH 08/11] Fix GMRES for complex operators 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. --- psydac/linalg/solvers.py | 17 +++++++++-------- psydac/linalg/tests/test_solvers.py | 6 ++++-- 2 files changed, 13 insertions(+), 10 deletions(-) diff --git a/psydac/linalg/solvers.py b/psydac/linalg/solvers.py index 82a459a5b..68cd732de 100644 --- a/psydac/linalg/solvers.py +++ b/psydac/linalg/solvers.py @@ -1791,7 +1791,7 @@ def solve(self, b, out=None): # update the residual vector beta.append(- sn[k] * beta[k]) - beta[k] *= cn[k] + beta[k] *= cn[k].conjugate() am = abs(beta[k+1]) if verbose: @@ -1834,7 +1834,7 @@ def arnoldi(self, k, p): self._A.dot( self._Q[k] , out=p) # Krylov vector for i in range(k + 1): # Modified Gram-Schmidt, keeping Hessenberg matrix - h[i] = p.inner(self._Q[i]) + h[i] = self._Q[i].inner(p) p.mul_iadd(-h[i], self._Q[i]) h[k+1] = sqrt(p.inner(p).real) @@ -1847,23 +1847,24 @@ def arnoldi(self, k, p): def apply_givens_rotation(self, k, sn, cn): # Apply Givens rotation to last column of H + # Rotation [[conj(c), conj(s)], [-s, c]] is unitary also for complex c, s h = self._H[:k+2, k] for i in range(k): h_i_prev = h[i] - h[i] *= cn[i] - h[i] += sn[i] * h[i+1] + h[i] *= cn[i].conjugate() + h[i] += sn[i].conjugate() * h[i+1] h[i+1] *= cn[i] h[i+1] -= sn[i] * h_i_prev - - mod = (h[k]**2 + h[k+1]**2)**0.5 + + mod = (abs(h[k])**2 + abs(h[k+1])**2)**0.5 cn.append( h[k] / mod ) sn.append( h[k+1] / mod ) - h[k] *= cn[k] - h[k] += sn[k] * h[k+1] + h[k] *= cn[k].conjugate() + h[k] += sn[k].conjugate() * h[k+1] h[k+1] = 0. # becomes triangular def dot(self, b, out=None): diff --git a/psydac/linalg/tests/test_solvers.py b/psydac/linalg/tests/test_solvers.py index 8edf38abc..e587519ee 100644 --- a/psydac/linalg/tests/test_solvers.py +++ b/psydac/linalg/tests/test_solvers.py @@ -270,10 +270,12 @@ def test_solver_diagonal(solver: str, diagonal: float, expected: float) -> None: #=============================================================================== # GMRES converges in at most n iterations, and must report the true residual @pytest.mark.parametrize('maxiter', [1, 6, 12]) -def test_GMRES_solve(maxiter: int) -> None: +@pytest.mark.parametrize('dtype', [float, complex]) +def test_GMRES_solve(maxiter: int, dtype: type) -> None: n = 12 - _, A, xe = define_data(n, 1, [-7, -6, -1]) + diagonals = [-7-2j, -6-2j, -1-10j] if dtype == complex else [-7, -6, -1] + _, A, xe = define_data(n, 1, diagonals, dtype=dtype) b = A @ xe solver = inverse(A, 'GMRES', tol=1e-10, maxiter=maxiter) From 899997d3f0be2a016036a06c4f15891184824138 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Mon, 5 Oct 2026 19:29:00 +0200 Subject: [PATCH 09/11] Fix preconditioned BiCGSTAB for 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`. --- psydac/linalg/solvers.py | 8 ++++---- psydac/linalg/tests/test_solvers.py | 4 +--- 2 files changed, 5 insertions(+), 7 deletions(-) diff --git a/psydac/linalg/solvers.py b/psydac/linalg/solvers.py index 68cd732de..c06a00b31 100644 --- a/psydac/linalg/solvers.py +++ b/psydac/linalg/solvers.py @@ -906,10 +906,10 @@ def solve_with_pc(self, b, out=None): while res_sqr > tol_sqr and niter < maxiter: - # v = A @ pp, vp = PC @ v, alphap = rhop/(vp.rp0) + # v = A @ pp, vp = PC @ v, alphap = rhop/(rp0.vp) A.dot(pp, out=v) pc.dot(v, out=vp) - alphap = rhop / vp.inner(rp0) + alphap = rhop / rp0.inner(vp) # s = r - alphap*v, sp = PC @ s r.copy(out=s) @@ -942,8 +942,8 @@ def solve_with_pc(self, b, out=None): tp *= omegap rp -= tp - # rhop_new = rp.rp0, betap = (alphap*rhop_new)/(omegap*rhop) - rhop_new = rp.inner(rp0) + # rhop_new = rp0.rp, betap = (alphap*rhop_new)/(omegap*rhop) + rhop_new = rp0.inner(rp) betap = (alphap*rhop_new) / (omegap*rhop) rhop = 1*rhop_new diff --git a/psydac/linalg/tests/test_solvers.py b/psydac/linalg/tests/test_solvers.py index e587519ee..f5e171710 100644 --- a/psydac/linalg/tests/test_solvers.py +++ b/psydac/linalg/tests/test_solvers.py @@ -85,9 +85,7 @@ def define_data(n, p, matrix_data, dtype=float): def test_solver_tridiagonal(n, p, dtype, solver, use_jacobi_pc, verbose=False): # Quickly skip tests that are not relevant - if solver == 'BiCGSTAB' and use_jacobi_pc and dtype == complex: - pytest.skip("Preconditioned BiCGSTAB only works for real matrices") - elif solver == 'MINRES' and dtype == complex: + if solver == 'MINRES' and dtype == complex: pytest.skip("MINRES only works for real matrices") # Also skip some problematic tests for now -- see Issue #557 From 8e829a1682b41b00a2d90098938417eaa80b7cc3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Mon, 5 Oct 2026 19:33:33 +0200 Subject: [PATCH 10/11] Fix CG crash with maxiter=1 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`. --- psydac/linalg/solvers.py | 6 ++++-- psydac/linalg/tests/test_solvers.py | 16 ++++++++++++++++ 2 files changed, 20 insertions(+), 2 deletions(-) diff --git a/psydac/linalg/solvers.py b/psydac/linalg/solvers.py index c06a00b31..443d61601 100644 --- a/psydac/linalg/solvers.py +++ b/psydac/linalg/solvers.py @@ -239,7 +239,8 @@ def solve_without_pc(self, b, out=None): template = "| {:7d} | {:19.2e} |" print(template.format(1, sqrt(am))) - # Iterate to convergence + # Iterate to convergence (iteration 1 is the initial residual) + m = 1 for m in range(2, maxiter+1): if am < tol_sqr: m -= 1 @@ -339,7 +340,8 @@ def solve_with_pc(self, b, out=None): template = "| {:7d} | {:19.2e} |" print( template.format(1, sqrt(nrmr_sqr))) - # Iterate to convergence + # Iterate to convergence (iteration 1 is the initial residual) + k = 1 for k in range(2, maxiter+1): if nrmr_sqr < tol_sqr: diff --git a/psydac/linalg/tests/test_solvers.py b/psydac/linalg/tests/test_solvers.py index f5e171710..b15a8cfb5 100644 --- a/psydac/linalg/tests/test_solvers.py +++ b/psydac/linalg/tests/test_solvers.py @@ -284,6 +284,22 @@ def test_GMRES_solve(maxiter: int, dtype: type) -> None: assert np.isclose(info['res_norm'], np.sqrt(r.inner(r).real), rtol=1e-8, atol=1e-12) assert info['success'] == (maxiter == n) +#=============================================================================== +# With maxiter=1, CG only evaluates the residual of the initial guess +@pytest.mark.parametrize('use_jacobi_pc', [False, True]) +def test_ConjugateGradient_solve_maxiter_1(use_jacobi_pc: bool) -> None: + + V, A, xe = define_data_hermitian(6, 1) + pc = A.diagonal(inverse=True) if use_jacobi_pc else None + + solver = inverse(A, 'CG', pc=pc, tol=1e-10, maxiter=1) + x = solver @ (A @ xe) + info = solver.get_info() + + assert np.array_equal(x.toarray(), V.zeros().toarray()) + assert info['niter'] == 1 + assert not info['success'] + #=============================================================================== def test_LST_preconditioner(comm=None): From caf0703d46a5884d367bf4c279429137a5d18143 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Mon, 5 Oct 2026 19:37:35 +0200 Subject: [PATCH 11/11] Fix LSMR division by zero when b = 0 or A^H r = 0 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`. --- psydac/linalg/solvers.py | 12 +++++++++++- psydac/linalg/tests/test_solvers.py | 18 ++++++++++++++++++ 2 files changed, 29 insertions(+), 1 deletion(-) diff --git a/psydac/linalg/solvers.py b/psydac/linalg/solvers.py index 443d61601..d0e6746b2 100644 --- a/psydac/linalg/solvers.py +++ b/psydac/linalg/solvers.py @@ -1492,6 +1492,15 @@ def solve(self, b, out=None): if conlim > 0:ctol = 1 / conlim normr = beta + # Early exit as in SciPy: the solution of A x = 0 is x = 0, and + # A^H (b - A x) = 0 means that x is already a least-squares solution + if normb == 0: + x *= 0.0 + normr = 0.0 + istop = 1 + elif alpha * beta == 0: + istop = 1 if beta == 0 else 2 + # Reverse the order here from the original matlab code because if verbose: @@ -1502,7 +1511,8 @@ def solve(self, b, out=None): template = "| {:7d} | {:19.2e} |" # Main iteration loop. - for itn in range(1, maxiter + 1): + while istop == 0 and itn < maxiter: + itn += 1 # Perform the next step of the bidiagonalization to obtain the # next beta, u, alpha, v. These satisfy the relations diff --git a/psydac/linalg/tests/test_solvers.py b/psydac/linalg/tests/test_solvers.py index b15a8cfb5..1eeab1fa4 100644 --- a/psydac/linalg/tests/test_solvers.py +++ b/psydac/linalg/tests/test_solvers.py @@ -300,6 +300,24 @@ def test_ConjugateGradient_solve_maxiter_1(use_jacobi_pc: bool) -> None: assert info['niter'] == 1 assert not info['success'] +#=============================================================================== +# LSMR returns before iterating if b = 0 (x = 0) or if A^H (b - A x0) = 0 (x = x0) +@pytest.mark.parametrize(('diagonal', 'rhs'), [(2.0, 0.0), (0.0, 1.0)]) +def test_LSMR_solve_early_exit(diagonal: float, rhs: float) -> None: + + V, A, xe = define_data(6, 1, [0.0, diagonal, 0.0]) + b = V.zeros() + b[:] = rhs + + solver = inverse(A, 'LSMR', x0=xe, tol=1e-10) + x = solver @ b + info = solver.get_info() + + expected = xe if rhs else V.zeros() + assert np.array_equal(x.toarray(), expected.toarray()) + assert info['niter'] == 0 + assert info['success'] + #=============================================================================== def test_LST_preconditioner(comm=None):