diff --git a/AGENTS.md b/AGENTS.md index 3390dfeed..84c2bc1ec 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -94,7 +94,9 @@ memory supercomputers using MPI and OpenMP. . Run `python -m isort --check $(git diff --name-only --diff-filter=d devel -- "*.py")` . Run `pylint --disable=all --enable=unused-import $(git diff --name-only --diff-filter=d devel -- "*.py")` - Until the code base is reformatted, black and isort fail on most existing files: - do not reformat code outside your change. The pylint check must pass. + do not reformat existing code, and do not run black on the changed lines only + (it mixes quote styles, comment separators, and alignment within a file). + The pylint check must pass. - 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) diff --git a/CHANGELOG.md b/CHANGELOG.md index cb92b760e..779ab50ae 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -25,6 +25,7 @@ All notable changes to this project will be documented in this file. ### Fixed +- #604 : Fix several bugs in iterative linear solvers (MINRES, LSMR, GMRES, CG, BiCGSTAB), see #603 - #527 : Require `sympde==0.20.0` which changes how multipatch interfaces are defined - #527 : Fix `psydac.cad.elevate` and `refine` for NURBS mappings, see #608 - #527 : Keep the NURBS weights in `psydac.cad.translate`, see #609 diff --git a/psydac/linalg/solvers.py b/psydac/linalg/solvers.py index fc3cce3cc..f966b31f1 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 @@ -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,12 @@ 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^H vp) + # `inner` conjugates its first argument, as `vdot` in SciPy's bicgstab: + # the arguments must be swapped if this convention changes (see #485) 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 +946,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^H rp, betap = (alphap*rhop_new)/(omegap*rhop) + rhop_new = rp0.inner(rp) betap = (alphap*rhop_new) / (omegap*rhop) rhop = 1*rhop_new @@ -1218,12 +1222,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 +1494,16 @@ def solve(self, b, out=None): if conlim > 0:ctol = 1 / conlim normr = beta + # Early exits as in SciPy (which returns istop = 0 in both cases): + # 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 +1514,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 +1616,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 +1797,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) @@ -1786,24 +1804,27 @@ def solve(self, b, out=None): # make the last diagonal entry in H equal to 0, so that H becomes upper triangular self.apply_givens_rotation(k, sn, cn) - # update the residual vector + # update the residual vector with the same rotation, as in PETSc 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 +1849,9 @@ 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] = Q[i]^H p, because `inner` conjugates its first argument: + # the arguments must be swapped if this convention changes (see #485) + h[i] = self._Q[i].inner(p) p.mul_iadd(-h[i], self._Q[i]) h[k+1] = sqrt(p.inner(p).real) @@ -1840,23 +1863,24 @@ def arnoldi(self, k, p): self._Q.append(p.copy()) def apply_givens_rotation(self, k, sn, cn): - # Apply Givens rotation to last column of H + # Apply Givens rotation [[conj(c), s], [-s, c]] to last column of H, + # as in PETSc: s is real, because h[k+1] is a norm (see `arnoldi`) h = self._H[:k+2, k] for i in range(k): h_i_prev = h[i] - h[i] *= cn[i] + h[i] *= cn[i].conjugate() h[i] += sn[i] * 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] *= cn[k].conjugate() h[k] += sn[k] * h[k+1] h[k+1] = 0. # becomes triangular diff --git a/psydac/linalg/tests/test_solvers.py b/psydac/linalg/tests/test_solvers.py index ee122d876..8683f1174 100644 --- a/psydac/linalg/tests/test_solvers.py +++ b/psydac/linalg/tests/test_solvers.py @@ -8,7 +8,7 @@ import pytest from mpi4py import MPI -from sympde.topology import elements_of, Mapping, Derham, Square, Cube +from sympde.topology import elements_of, Derham, Square, Cube from sympde.calculus import inner from sympde.expr import integral, BilinearForm @@ -24,31 +24,42 @@ def define_data_hermitian(n, p, dtype=float): - domain_decomposition = DomainDecomposition([n - p], [False]) - cart = CartDecomposition(domain_decomposition, [n], [np.array([0])], [np.array([n - 1])], [p], [1]) - # ... Vector Spaces - V = StencilVectorSpace(cart,dtype=dtype) - e = V.ends[0] - 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 - if dtype==complex: - factor=1+1j - else: - factor=1 - A = StencilMatrix(V, V) - A[:, -p:0] = 1-1*factor - A[:, 0:1] = 2 * p - A[:, 1:p + 1] = 1-1*factor.conjugate() - A.remove_spurious_entries() - - # Build exact solution - xe = StencilVector(V) - xe[s:e + 1] = factor*np.random.random(e + 1 - s) - return(V, A, xe) + """ + Create a Hermitian positive definite banded matrix, and a random exact solution. + + The matrix is created by `define_data` with 2p on the main diagonal, and -1 + (real case) or -1j (complex case) on the p lower diagonals. The p upper + diagonals hold the complex conjugate of the lower ones. In the complex case, + the exact solution is complex as well. See `define_data` for the returns. + """ + lower = -1j if dtype == complex else -1 + V, A, xe = define_data(n, p, [lower, 2 * p, np.conj(lower)], dtype=dtype) + if dtype == complex: + xe *= 1 + 1j + return V, A, xe def define_data(n, p, matrix_data, dtype=float): + """ + Create a banded matrix with 2p+1 constant diagonals, and a random exact solution. + + The p lower diagonals are set to matrix_data[0], the main diagonal to + matrix_data[1], and the p upper diagonals to matrix_data[2]. Hence the + matrix is not necessarily symmetric or positive definite. It is Hermitian + positive definite if matrix_data[2] = conj(matrix_data[0]) and + matrix_data[1] >= 2p |matrix_data[0]|, e.g. with matrix_data = [-1, 2p, -1]. + See also `define_data_hermitian`. + + Returns + ------- + V : StencilVectorSpace + The space of the matrix and of the solution. + + A : StencilMatrix + The banded matrix. + + xe : StencilVector + The exact solution, with random entries. + """ domain_decomposition = DomainDecomposition([n - p], [False]) cart = CartDecomposition(domain_decomposition, [n], [np.array([0])], [np.array([n - 1])], [p], [1]) # ... Vector Spaces @@ -56,11 +67,8 @@ def define_data(n, p, matrix_data, dtype=float): e = V.ends[0] 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 - A = StencilMatrix(V, V) - A[:, -p:0] = -matrix_data[0] + A[:, -p:0] = matrix_data[0] A[:, 0:1] = matrix_data[1] A[:, 1:p + 1] = matrix_data[2] A.remove_spurious_entries() @@ -85,9 +93,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 @@ -107,10 +113,10 @@ def test_solver_tridiagonal(n, p, dtype, solver, use_jacobi_pc, verbose=False): diagonals = None elif solver in ['BiCG', 'BiCGSTAB', 'LSMR']: hermitian = False - diagonals = [1-10j, 6+9j, 3+5j] if dtype==complex else [1, 6, 3] + diagonals = [-1+10j, 6+9j, 3+5j] if dtype==complex else [-1, 6, 3] elif solver == 'GMRES': hermitian = False - diagonals = [-7-2j, -6-2j, -1-10j] if dtype==complex else [-7, -6, -1] + diagonals = [7+2j, -6-2j, -1-10j] if dtype==complex else [7, -6, -1] if hermitian: V, A, xe = define_data_hermitian(n, p, dtype=dtype) @@ -252,6 +258,82 @@ 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 +#=============================================================================== +# With A = 0, the stopping tests of MINRES set test1 = test2 = inf +def test_MinimumResidual_solve_zero_operator() -> None: + + V, A, _ = define_data(6, 1, [0.0, 0.0, 0.0]) + b = V.zeros() + b[:] = 1.0 + + x = inverse(A, 'MINRES', tol=1e-10) @ b + + assert np.array_equal(x.toarray(), np.zeros(6)) + +#=============================================================================== +# With A = 2I, LSMR finds the exact solution, and its stopping tests set test2 = inf +def test_LSMR_solve_zero_residual() -> None: + + V, A, _ = define_data(6, 1, [0.0, 2.0, 0.0]) + b = V.zeros() + b[:] = 1.0 + + x = inverse(A, 'LSMR', tol=1e-10) @ b + + assert np.array_equal(x.toarray(), np.full(6, 0.5)) + +#=============================================================================== +# GMRES must report the true residual of the returned solution, converged or not +@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) + +#=============================================================================== +# 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):