Fix several bugs in iterative linear solvers

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
 

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 `<module>_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_<module>.py`
+- Test function names follow `test_<function>` or `test_<class>_<method>`
+- 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`

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

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
 #==============================================================================

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
@@ -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:
@@ -906,10 +908,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 +944,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
 
@@ -1218,12 +1220,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
@@ -1485,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:
@@ -1495,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
@@ -1596,7 +1613,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
@@ -1777,8 +1794,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)
@@ -1788,22 +1803,25 @@ 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:
                 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"])
@@ -1828,7 +1846,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)
@@ -1841,23 +1859,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):

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
@@ -252,6 +250,74 @@ 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))
+
+#===============================================================================
+# GMRES converges in at most n iterations, and must report the true residual
+@pytest.mark.parametrize('maxiter', [1, 6, 12])
+@pytest.mark.parametrize('dtype', [float, complex])
+def test_GMRES_solve(maxiter: int, dtype: type) -> None:
+
+    n = 12
+    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)
+    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)
+
+#===============================================================================
+# 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']
+
+#===============================================================================
+# 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):