diff --git a/.github/workflows/continuous-integration.yml b/.github/workflows/continuous-integration.yml index 826027fa1..95c074f5f 100644 --- a/.github/workflows/continuous-integration.yml +++ b/.github/workflows/continuous-integration.yml @@ -16,11 +16,10 @@ jobs: fail-fast: false matrix: os: [ ubuntu-latest ] - python-version: [ 3.8, 3.9, '3.10', '3.11' ] + python-version: [ 3.9, '3.10', '3.11' ] isMerge: - ${{ github.event_name == 'push' && github.ref == 'refs/heads/devel' }} exclude: - - { isMerge: false, python-version: 3.9 } - { isMerge: false, python-version: '3.10' } include: - os: macos-latest diff --git a/psydac/api/settings.py b/psydac/api/settings.py index 6d0988451..74a48c645 100644 --- a/psydac/api/settings.py +++ b/psydac/api/settings.py @@ -67,4 +67,4 @@ 'pyccel-intel' : PSYDAC_BACKEND_IPYCCEL, 'pyccel-pgi' : PSYDAC_BACKEND_PGPYCCEL, 'pyccel-nvidia': PSYDAC_BACKEND_NVPYCCEL, -} \ No newline at end of file +} diff --git a/psydac/api/tests/test_api_2d_compatible_spaces.py b/psydac/api/tests/test_api_2d_compatible_spaces.py index feb7cd439..45c04a943 100644 --- a/psydac/api/tests/test_api_2d_compatible_spaces.py +++ b/psydac/api/tests/test_api_2d_compatible_spaces.py @@ -130,7 +130,7 @@ def run_stokes_2d_dir(domain, f, ue, pe, *, homogeneous, ncells, degree, scipy=F # ... solve linear system using scipy.sparse.linalg or psydac if scipy: - tol = 1e-11 + rtol = 1e-11 equation_h.assemble() A0 = equation_h.linear_system.lhs.tosparse() b0 = equation_h.linear_system.rhs.toarray() @@ -145,17 +145,17 @@ def run_stokes_2d_dir(domain, f, ue, pe, *, homogeneous, ncells, degree, scipy=F A1 = a1_h.assemble().tosparse() b1 = l1_h.assemble().toarray() - x1, info = sp_minres(A1, b1, tol=tol) + x1, info = sp_minres(A1, b1, rtol=rtol) print('Boundary solution with scipy.sparse: success = {}'.format(info == 0)) - x0, info = sp_minres(A0, b0 - A0.dot(x1), tol=tol) + x0, info = sp_minres(A0, b0 - A0.dot(x1), rtol=rtol) print('Interior solution with scipy.sparse: success = {}'.format(info == 0)) # Solution is sum of boundary and interior contributions x = x0 + x1 else: - x, info = sp_minres(A0, b0, tol=tol) + x, info = sp_minres(A0, b0, rtol=rtol) print('Solution with scipy.sparse: success = {}'.format(info == 0)) # Convert to stencil format diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py index 81079a5c7..33fba4967 100644 --- a/psydac/linalg/basic.py +++ b/psydac/linalg/basic.py @@ -8,6 +8,8 @@ """ from abc import ABC, abstractmethod +from types import LambdaType +from inspect import signature import numpy as np from scipy.sparse import coo_matrix @@ -25,7 +27,8 @@ 'ComposedLinearOperator', 'PowerLinearOperator', 'InverseLinearOperator', - 'LinearSolver' + 'LinearSolver', + 'MatrixFreeLinearOperator' ) #=============================================================================== @@ -1111,3 +1114,108 @@ def solve(self, rhs, out=None): @property def T(self): return self.transpose() + +#=============================================================================== +class MatrixFreeLinearOperator(LinearOperator): + """ + General linear operator represented by a callable dot method. + + Parameters + ---------- + domain : VectorSpace + The domain of the linear operator. + + codomain : VectorSpace + The codomain of the linear operator. + + dot : Callable + The method of the linear operator, assumed to map from domain to codomain. + This method can take out as an optional argument but this is not mandatory. + + dot_transpose: Callable + The method of the transpose of the linear operator, assumed to map from codomain to domain. + This method can take out as an optional argument but this is not mandatory. + + Examples + -------- + # example 1: a matrix encapsulated as a (fake) matrix-free linear operator + A_SM = StencilMatrix(V, W) + AT_SM = A_SM.transpose() + A = MatrixFreeLinearOperator(domain=V, codomain=W, dot=lambda v: A_SM @ v, dot_transpose=lambda v: AT_SM @ v) + + # example 2: a truly matrix-free linear operator + A = MatrixFreeLinearOperator(domain=V, codomain=V, dot=lambda v: 2*v, dot_transpose=lambda v: 2*v) + + """ + + def __init__(self, domain, codomain, dot, dot_transpose=None): + + assert isinstance(domain, VectorSpace) + assert isinstance(codomain, VectorSpace) + assert isinstance(dot, LambdaType) + + self._domain = domain + self._codomain = codomain + self._dot = dot + + sig = signature(dot) + self._dot_takes_out_arg = ('out' in [p.name for p in sig.parameters.values() if p.kind == p.KEYWORD_ONLY]) + + if dot_transpose is not None: + assert isinstance(dot_transpose, LambdaType) + self._dot_transpose = dot_transpose + sig = signature(dot_transpose) + self._dot_transpose_takes_out_arg = ('out' in [p.name for p in sig.parameters.values() if p.kind == p.KEYWORD_ONLY]) + else: + self._dot_transpose = None + self._dot_transpose_takes_out_arg = False + + @property + def domain(self): + return self._domain + + @property + def codomain(self): + return self._codomain + + @property + def dtype(self): + return None + + def dot(self, v, out=None): + assert isinstance(v, Vector) + assert v.space == self.domain + + if out is not None: + assert isinstance(out, Vector) + assert out.space == self.codomain + else: + out = self.codomain.zeros() + + if self._dot_takes_out_arg: + self._dot(v, out=out) + else: + # provided dot product does not take an out argument: we simply copy the result into out + self._dot(v).copy(out=out) + + return out + + def toarray(self): + raise NotImplementedError('toarray() is not defined for MatrixFreeLinearOperator.') + + def tosparse(self): + raise NotImplementedError('tosparse() is not defined for MatrixFreeLinearOperator.') + + def transpose(self, conjugate=False): + if self._dot_transpose is None: + raise NotImplementedError('no transpose dot method was given -- cannot create the transpose operator') + + if conjugate: + if self._dot_transpose_takes_out_arg: + new_dot = lambda v, out=None: self._dot_transpose(v, out=out).conjugate() + else: + new_dot = lambda v: self._dot_transpose(v).conjugate() + else: + new_dot = self._dot_transpose + + return MatrixFreeLinearOperator(domain=self.codomain, codomain=self.domain, dot=new_dot, dot_transpose=self._dot) diff --git a/psydac/linalg/solvers.py b/psydac/linalg/solvers.py index 297720ab7..862868ea2 100644 --- a/psydac/linalg/solvers.py +++ b/psydac/linalg/solvers.py @@ -39,7 +39,7 @@ def inverse(A, solver, **kwargs): A : psydac.linalg.basic.LinearOperator Left-hand-side matrix A of linear system; individual entries A[i,j] can't be accessed, but A has 'shape' attribute and provides 'dot(p)' - function (i.e. matrix-vector product A*p). + function (e.g. a matrix-vector product A*p). solver : str Preferred iterative solver. Options are: 'cg', 'pcg', 'bicg', @@ -1165,7 +1165,7 @@ def solve(self, b, out=None): A.dot(x, out=y) y -= b y *= -1.0 - y.copy(out=res_old) + y.copy(out=res_old) # res = b - A*x beta = sqrt(res_old.dot(res_old)) @@ -1193,8 +1193,15 @@ def solve(self, b, out=None): print( "+---------+---------------------+") template = "| {:7d} | {:19.2e} |" + # check whether solution is already converged: + if beta < tol: + istop = 1 + rnorm = beta + if verbose: + print( template.format(itn, rnorm )) - for itn in range(1, maxiter + 1 ): + while istop == 0 and itn < maxiter: + itn += 1 s = 1.0/beta y.copy(out=v) diff --git a/psydac/linalg/tests/test_linalg.py b/psydac/linalg/tests/test_linalg.py index 8b30802b8..15d5ad312 100644 --- a/psydac/linalg/tests/test_linalg.py +++ b/psydac/linalg/tests/test_linalg.py @@ -20,7 +20,7 @@ def array_equal(a, b): def sparse_equal(a, b): return (a.tosparse() != b.tosparse()).nnz == 0 -def is_pos_def(A): +def assert_pos_def(A): assert isinstance(A, LinearOperator) A_array = A.toarray() assert np.all(np.linalg.eigvals(A_array) > 0) @@ -50,7 +50,7 @@ def get_StencilVectorSpace(n1, n2, p1, p2, P1, P2): V = StencilVectorSpace(C) return V -def get_positive_definite_stencilmatrix(V): +def get_positive_definite_StencilMatrix(V): np.random.seed(2) assert isinstance(V, StencilVectorSpace) @@ -700,9 +700,9 @@ def test_positive_definite_matrix(n1, n2, p1, p2): P1 = False P2 = False V = get_StencilVectorSpace(n1, n2, p1, p2, P1, P2) - S = get_positive_definite_stencilmatrix(V) + S = get_positive_definite_StencilMatrix(V) - is_pos_def(S) + assert_pos_def(S) #=============================================================================== @pytest.mark.parametrize('n1', [3, 5]) @@ -745,7 +745,7 @@ def test_operator_evaluation(n1, n2, p1, p2): V = get_StencilVectorSpace(n1, n2, p1, p2, P1, P2) # Initiate positive definite StencilMatrices for which the cg inverse works (necessary for certain tests) - S = get_positive_definite_stencilmatrix(V) + S = get_positive_definite_StencilMatrix(V) # Initiate StencilVectors v = StencilVector(V) @@ -769,7 +769,7 @@ def test_operator_evaluation(n1, n2, p1, p2): ### 2.1 PowerLO test Bmat = B.toarray() - is_pos_def(B) + assert_pos_def(B) uarr = u.toarray() b0 = ( B**0 @ u ).toarray() b1 = ( B**1 @ u ).toarray() @@ -799,7 +799,7 @@ def test_operator_evaluation(n1, n2, p1, p2): assert np.array_equal(zeros, z2) Smat = S.toarray() - is_pos_def(S) + assert_pos_def(S) varr = v.toarray() s0 = ( S**0 @ v ).toarray() s1 = ( S**1 @ v ).toarray() @@ -960,8 +960,8 @@ def test_x0update(solver): P1 = False P2 = False V = get_StencilVectorSpace(n1, n2, p1, p2, P1, P2) - A = get_positive_definite_stencilmatrix(V) - is_pos_def(A) + A = get_positive_definite_StencilMatrix(V) + assert_pos_def(A) b = StencilVector(V) for n in range(n1): b[n, :] = 1. diff --git a/psydac/linalg/tests/test_matrix_free.py b/psydac/linalg/tests/test_matrix_free.py new file mode 100644 index 000000000..e475fb15f --- /dev/null +++ b/psydac/linalg/tests/test_matrix_free.py @@ -0,0 +1,127 @@ +import pytest +import numpy as np + +from psydac.linalg.block import BlockLinearOperator, BlockVector, BlockVectorSpace +from psydac.linalg.basic import LinearOperator, ZeroOperator, IdentityOperator, ComposedLinearOperator, SumLinearOperator, PowerLinearOperator, ScaledLinearOperator +from psydac.linalg.basic import MatrixFreeLinearOperator +from psydac.linalg.stencil import StencilVectorSpace, StencilVector, StencilMatrix +from psydac.linalg.solvers import ConjugateGradient, inverse +from psydac.ddm.cart import DomainDecomposition, CartDecomposition + +from psydac.linalg.tests.test_linalg import get_StencilVectorSpace, get_positive_definite_StencilMatrix, assert_pos_def + +def get_random_StencilMatrix(domain, codomain): + + np.random.seed(2) + V = domain + W = codomain + assert isinstance(V, StencilVectorSpace) + assert isinstance(W, StencilVectorSpace) + [n1, n2] = V._npts + [p1, p2] = V._pads + [P1, P2] = V._periods + assert (P1 == False) and (P2 == False) + + [m1, m2] = W._npts + [q1, q2] = W._pads + [Q1, Q2] = W._periods + assert (Q1 == False) and (Q2 == False) + + S = StencilMatrix(V, W) + + for i in range(0, q1+1): + if i != 0: + for j in range(-q2, q2+1): + S[:, :, i, j] = 2*np.random.random()-1 + else: + for j in range(1, q2+1): + S[:, :, i, j] = 2*np.random.random()-1 + S.remove_spurious_entries() + + return S + +def get_random_StencilVector(V): + np.random.seed(3) + assert isinstance(V, StencilVectorSpace) + [n1, n2] = V._npts + v = StencilVector(V) + for i in range(n1): + for j in range(n2): + v[i,j] = np.random.random() + return v + +#=============================================================================== +@pytest.mark.parametrize('n1', [3, 5]) +@pytest.mark.parametrize('n2', [4, 7]) +@pytest.mark.parametrize('p1', [2, 6]) +@pytest.mark.parametrize('p2', [3, 9]) + +def test_fake_matrix_free(n1, n2, p1, p2): + P1 = False + P2 = False + m1 = (n2+n1)//2 + m2 = n1+1 + q1 = p1 # using same degrees because both spaces must have same padding for now + q2 = p2 + V1 = get_StencilVectorSpace(n1, n2, p1, p2, P1, P2) + V2 = get_StencilVectorSpace(m1, m2, q1, q2, P1, P2) + S = get_random_StencilMatrix(codomain=V2, domain=V1) + O = MatrixFreeLinearOperator(codomain=V2, domain=V1, dot=lambda v: S @ v) + + print(f'O.domain = {O.domain}') + print(f'S.domain = {S.domain}') + print(f'V1: = {V1}') + v = get_random_StencilVector(V1) + tol = 1e-10 + y = S.dot(v) + x = O.dot(v) + print(f'error = {np.linalg.norm( (x - y).toarray() )}') + assert np.linalg.norm( (x - y).toarray() ) < tol + O.dot(v, out=x) + print(f'error = {np.linalg.norm( (x - y).toarray() )}') + assert np.linalg.norm( (x - y).toarray() ) < tol + +@pytest.mark.parametrize('solver', ['cg', 'pcg', 'bicg', 'minres', 'lsmr']) + +def test_solvers_matrix_free(solver): + print(f'solver = {solver}') + n1 = 4 + n2 = 3 + p1 = 5 + p2 = 2 + P1 = False + P2 = False + V = get_StencilVectorSpace(n1, n2, p1, p2, P1, P2) + A_SM = get_positive_definite_StencilMatrix(V) + assert_pos_def(A_SM) + AT_SM = A_SM.transpose() + A = MatrixFreeLinearOperator(domain=V, codomain=V, dot=lambda v: A_SM @ v, dot_transpose=lambda v: AT_SM @ v) + + # get rhs and solution + b = get_random_StencilVector(V) + x = A.dot(b) + + # Create Inverse with A + tol = 1e-6 + if solver == 'pcg': + inv_diagonal = A_SM.diagonal(inverse=True) + A_inv = inverse(A, solver, pc=inv_diagonal, tol=tol) + else: + A_inv = inverse(A, solver, tol=tol) + + AA = A_inv._A + xx = AA.dot(b) + print(f'norm(xx) = {np.linalg.norm( xx.toarray() )}') + print(f'norm(x) = {np.linalg.norm( x.toarray() )}') + + # Apply inverse and check + y = A_inv @ x + error = np.linalg.norm( (b - y).toarray()) + assert np.linalg.norm( (b - y).toarray() ) < tol + +#=============================================================================== +# SCRIPT FUNCTIONALITY +#=============================================================================== +if __name__ == "__main__": + import sys + pytest.main( sys.argv ) diff --git a/pyproject.toml b/pyproject.toml index d59494d00..a797463b5 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -7,7 +7,7 @@ name = "psydac" version = "0.1" description = "Python package for isogeometric analysis (IGA)" readme = "README.md" -requires-python = ">= 3.8, < 3.12" +requires-python = ">= 3.9, < 3.12" license = {file = "LICENSE"} authors = [ {name = "Psydac development team", email = "psydac@googlegroups.com"} @@ -23,7 +23,7 @@ dependencies = [ # Third-party packages from PyPi 'numpy >= 1.16', - 'scipy >= 0.18, < 1.14', + 'scipy >= 1.12', 'sympy >= 1.5', 'matplotlib', 'pytest >= 4.5', @@ -50,10 +50,7 @@ dependencies = [ 'tblib', # IGAKIT - not on PyPI - - # !! WARNING !! Path to igakit below is from fork pyccel/igakit. This was done to - # quickly fix the numpy 2.0 issue. See https://github.com/dalcinl/igakit/pull/4 - 'igakit @ https://github.com/pyccel/igakit/archive/refs/heads/bugfix-numpy2.0.zip' + 'igakit @ https://github.com/dalcinl/igakit/archive/refs/heads/master.zip' ] [project.urls] diff --git a/requirements.txt b/requirements.txt index 95b01680d..a42e47520 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,6 +1,7 @@ wheel setuptools >= 61, != 67.2.0 numpy >= 1.16 +scipy >= 1.12 Cython >= 0.25, < 3.0 mpi4py >= 3.0.0, < 4