diff --git a/.github/workflows/continuous-integration.yml b/.github/workflows/continuous-integration.yml index ae4c2e09d..f9f7b9636 100644 --- a/.github/workflows/continuous-integration.yml +++ b/.github/workflows/continuous-integration.yml @@ -179,13 +179,13 @@ jobs: --cov psydac --cov-config $GITHUB_WORKSPACE/pyproject.toml --cov-report xml - --pyargs psydac -m "not parallel and not petsc" + --pyargs psydac.linalg -m "not parallel and not petsc" - name: Run single-process tests with Pytest on Ubuntu if: matrix.os == 'ubuntu-24.04' working-directory: ./pytest run: | - python -m pytest -n auto --pyargs psydac -m "not parallel and not petsc" + python -m pytest -n auto --pyargs psydac.linalg -m "not parallel and not petsc" - name: Upload coverage report to Codacy if: matrix.os == 'macos-14' @@ -203,17 +203,17 @@ jobs: - name: Run MPI tests with Pytest working-directory: ./pytest run: | - python mpi_tester.py --mpirun="mpiexec -n 4 ${MPI_OPTS}" --pyargs psydac -m "parallel and not petsc" + python mpi_tester.py --mpirun="mpiexec -n 4 ${MPI_OPTS}" --pyargs psydac.linalg -m "parallel and not petsc" - name: Run single-process PETSc tests with Pytest working-directory: ./pytest run: | - python -m pytest -n auto --pyargs psydac -m "not parallel and petsc" + python -m pytest -n auto --pyargs psydac.linalg -m "not parallel and petsc" - name: Run MPI PETSc tests with Pytest working-directory: ./pytest run: | - python mpi_tester.py --mpirun="mpiexec -n 4 ${MPI_OPTS}" --pyargs psydac -m "parallel and petsc" + python mpi_tester.py --mpirun="mpiexec -n 4 ${MPI_OPTS}" --pyargs psydac.linalg -m "parallel and petsc" - name: Remove test directory if: always() diff --git a/psydac/feec/multipatch/fem_linear_operators.py b/psydac/feec/multipatch/fem_linear_operators.py index aafb77931..7ffa8d9fe 100644 --- a/psydac/feec/multipatch/fem_linear_operators.py +++ b/psydac/feec/multipatch/fem_linear_operators.py @@ -58,6 +58,9 @@ def T(self): @property def dtype( self ): return self.domain.dtype + + def copy(self): + return FemLinearOperator(self.fem_domain, self.fem_codomain, self.matrix, self._sparse_matrix) def toarray(self): return self._matrix.toarray() @@ -135,6 +138,9 @@ def __init__( self, operators ): # matrix not defined by matrix product because it could break the Stencil Matrix structure + def copy(self): + return ComposedLinearOperator(self._operators) + def to_sparse_matrix( self, **kwargs): mat = self._operators[-1].to_sparse_matrix() for i in range(2, self._n+1): @@ -160,6 +166,9 @@ class IdLinearOperator( FemLinearOperator ): def __init__( self, V ): FemLinearOperator.__init__(self, fem_domain=V) + def copy(self): + return IdLinearOperator(self.fem_domain) + def to_sparse_matrix( self , **kwargs): return sparse_id( self.fem_domain.nbasis ) @@ -183,6 +192,9 @@ def __init__( self, B, A ): self._A = A self._B = B + def copy(self): + return SumLinearOperator(self._B, self._A) + def to_sparse_matrix( self, **kwargs): return self._A.to_sparse_matrix() + self._B.to_sparse_matrix() @@ -205,6 +217,9 @@ def __init__( self, c, A ): self._A = A self._c = c + def copy(self): + return MultLinearOperator(self._c, self._A) + def to_sparse_matrix( self, **kwargs): return self._c * self._A.to_sparse_matrix() diff --git a/psydac/feec/multipatch/operators.py b/psydac/feec/multipatch/operators.py index 0ee00210c..3fff465b0 100644 --- a/psydac/feec/multipatch/operators.py +++ b/psydac/feec/multipatch/operators.py @@ -472,6 +472,9 @@ def __init__( storage_fn) save_npz(storage_fn, self._sparse_matrix) + def copy(self): + raise NotImplementedError + def set_homogenous_bc(self, boundary, rhs=None): domain = self.symbolic_domain Vh = self.fem_domain @@ -721,6 +724,9 @@ def __init__( storage_fn) save_npz(storage_fn, self._sparse_matrix) + def copy(self): + raise NotImplementedError + def set_homogenous_bc(self, boundary): domain = self.symbolic_domain Vh = self.fem_domain @@ -1011,6 +1017,9 @@ def __init__( # matrices are not stored, we will probably compute them later pass + def copy(self): + raise NotImplementedError + def to_sparse_matrix(self): """ the Hodge matrix is the patch-wise multi-patch mass matrix @@ -1098,6 +1107,9 @@ def __init__(self, V0h, V1h): def transpose(self, conjugate=False): # todo (MCP): define as the dual differential operator return BrokenTransposedGradient_2D(self.fem_domain, self.fem_codomain) + + def copy(self): + return BrokenGradient_2D(self.fem_domain, self.fem_codomain) # ============================================================================== @@ -1116,6 +1128,9 @@ def __init__(self, V0h, V1h): def transpose(self, conjugate=False): # todo (MCP): discard return BrokenGradient_2D(self.fem_codomain, self.fem_domain) + + def copy(self): + return BrokenTransposedGradient_2D(self.fem_domain, self.fem_codomain) # ============================================================================== @@ -1132,6 +1147,9 @@ def __init__(self, V1h, V2h): def transpose(self, conjugate=False): return BrokenTransposedScalarCurl_2D( V1h=self.fem_domain, V2h=self.fem_codomain) + + def copy(self): + return BrokenScalarCurl_2D(self.fem_domain, self.fem_codomain) # ============================================================================== @@ -1148,6 +1166,9 @@ def __init__(self, V1h, V2h): def transpose(self, conjugate=False): return BrokenScalarCurl_2D(V1h=self.fem_codomain, V2h=self.fem_domain) + + def copy(self): + return BrokenTransposedScalarCurl_2D(self.fem_domain, self.fem_codomain) # ============================================================================== diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py index bf88dbd07..556b19ee2 100644 --- a/psydac/linalg/basic.py +++ b/psydac/linalg/basic.py @@ -280,6 +280,10 @@ def codomain(self): def dtype(self): pass + @abstractmethod + def copy(self): + pass + @abstractmethod def tosparse(self): pass @@ -303,8 +307,6 @@ def transpose(self, conjugate=False): """ pass - # TODO: check if we should add a copy method!!! - #------------------------------------- # Magic methods #------------------------------------- @@ -616,6 +618,14 @@ def operator(self): @property def dtype(self): return None + + def copy(self): + return ScaledLinearOperator(self.domain, self.codomain, self.scalar, self.operator.copy()) + + def set_scalar(self, c): + assert np.isscalar(c) + assert np.iscomplexobj(c) == (self.codomain._dtype == complex) + self._scalar = c def toarray(self): return self._scalar*self._operator.toarray() @@ -699,6 +709,10 @@ def addends(self): @property def dtype(self): return None + + def copy(self): + copied_addends = [addend.copy() for addend in self.addends] + return SumLinearOperator(self.domain, self.codomain, *copied_addends) def toarray(self): out = np.zeros(self.shape, dtype=self.dtype) @@ -830,6 +844,10 @@ def multiplicants(self): @property def dtype(self): return None + + def copy(self): + copied_multiplicants = [multiplicant.copy() for multiplicant in self.multiplicants] + return ComposedLinearOperator(self.domain, self.codomain, *copied_multiplicants) def toarray(self): raise NotImplementedError('toarray() is not defined for ComposedLinearOperators.') @@ -937,6 +955,9 @@ def operator(self): def factorial(self): """ Returns the power to which the operator is raised. """ return self._factorial + + def copy(self): + return PowerLinearOperator(self.domain, self.codomain, self.operator.copy(), self.factorial) def toarray(self): raise NotImplementedError('toarray() is not defined for PowerLinearOperators.') @@ -1201,6 +1222,9 @@ def codomain(self): @property def dtype(self): return None + + def copy(self): + return MatrixFreeLinearOperator(self.domain, self.codomain, self._dot, self._dot_transpose) def dot(self, v, out=None, **kwargs): assert isinstance(v, Vector) diff --git a/psydac/linalg/block.py b/psydac/linalg/block.py index f2a755015..7c5af0a7e 100644 --- a/psydac/linalg/block.py +++ b/psydac/linalg/block.py @@ -767,8 +767,8 @@ def __add__(self, M): for ij in set(self._blocks.keys()) | set(M._blocks.keys()): Bij = self[ij] Mij = M[ij] - if Bij is None: blocks[ij] = Mij.copy() - elif Mij is None: blocks[ij] = Bij.copy() + if Bij is None: blocks[ij] = Mij + elif Mij is None: blocks[ij] = Bij else : blocks[ij] = Bij + Mij mat = BlockLinearOperator(self.domain, self.codomain, blocks=blocks) if len(mat._blocks) != len(self._blocks): @@ -782,6 +782,7 @@ def __add__(self, M): # ... def __sub__(self, M): + from psydac.linalg.basic import ZeroOperator if not isinstance(M, BlockLinearOperator): return LinearOperator.__sub__(self, M) @@ -792,7 +793,7 @@ def __sub__(self, M): Bij = self[ij] Mij = M[ij] if Bij is None: blocks[ij] = -Mij - elif Mij is None: blocks[ij] = Bij.copy() + elif Mij is None: blocks[ij] = Bij else : blocks[ij] = Bij - Mij mat = BlockLinearOperator(self.domain, self.codomain, blocks=blocks) if len(mat._blocks) != len(self._blocks): @@ -1001,54 +1002,61 @@ def copy(self, out=None): return out # ... - def __imul__(self, a): - for Bij in self._blocks.values(): - Bij *= a - return self + # def __imul__(self, a): + # for Bij in self._blocks.values(): + # Bij *= a + # return self # ... - def __iadd__(self, M): - if not isinstance(M, BlockLinearOperator): - return LinearOperator.__add__(self, M) - - assert M. domain is self. domain - assert M.codomain is self.codomain - - for ij in set(self._blocks.keys()) | set(M._blocks.keys()): - - Mij = M[ij] - if Mij is None: - continue - - Bij = self[ij] - if Bij is None: - self[ij] = Mij.copy() - else: - Bij += Mij - - return self + # def __iadd__(self, M): + # if not isinstance(M, BlockLinearOperator): + # return LinearOperator.__add__(self, M) + + # assert M. domain is self. domain + # assert M.codomain is self.codomain + + # for ij in set(self._blocks.keys()) | set(M._blocks.keys()): + + # Mij = M[ij] + # if Mij is None: + # continue + + # Bij = self[ij] + # if Bij is None: + # Need to investigate: Removing .copy() here causes + # M = BlockLinearOperator(W, W, blocks=[[M1, M2], [M3, None]]) + # A = BlockLinearOperator(W, W) + # A += M + # A -= 2*M + # assert np.allclose(A.blocks[0][0]._data, -(M1)._data, rtol=1e-14, atol=1e-14 ) + # to fail. Possible explanation "-= 2*M" caues the entries of A to be multiplied as well, hence the result would be 0. + # self[ij] = Mij.copy() + # else: + # Bij += Mij + + # return self # ... - def __isub__(self, M): - if not isinstance(M, BlockLinearOperator): - return LinearOperator.__sub__(self, M) + # def __isub__(self, M): + # if not isinstance(M, BlockLinearOperator): + # return LinearOperator.__sub__(self, M) - assert M. domain is self. domain - assert M.codomain is self.codomain + # assert M. domain is self. domain + # assert M.codomain is self.codomain - for ij in set(self._blocks.keys()) | set(M._blocks.keys()): + # for ij in set(self._blocks.keys()) | set(M._blocks.keys()): - Mij = M[ij] - if Mij is None: - continue + # Mij = M[ij] + # if Mij is None: + # continue - Bij = self[ij] - if Bij is None: - self[ij] = -Mij - else: - Bij -= Mij + # Bij = self[ij] + # if Bij is None: + # self[ij] = -Mij + # else: + # Bij -= Mij - return self + # return self # ... def topetsc(self): @@ -1284,7 +1292,8 @@ def set_backend(self, backend): if not all(isinstance(b, (StencilMatrix, StencilInterfaceMatrix)) for b in self._blocks.values()): for b in self._blocks.values(): - b.set_backend(backend) + if getattr(b, "set_backend", None) is not None: + b.set_backend(backend) return block_shape = (self.n_block_rows, self.n_block_cols) diff --git a/psydac/linalg/fft.py b/psydac/linalg/fft.py index 5d67d178b..dbe87a695 100644 --- a/psydac/linalg/fft.py +++ b/psydac/linalg/fft.py @@ -121,9 +121,14 @@ def __init__(self, space, norm=None, workers=os.environ.get('OMP_NUM_THREADS', N assert isinstance(space, StencilVectorSpace) assert np.dtype(space.dtype).kind == 'c' workers = int(workers) if workers is not None else None + self._workers = workers + self._norm = norm super().__init__(space, lambda out: scifft.fft( out, axis=1, overwrite_x=True, workers=workers, norm=norm)) + + def copy(self): + return DistributedFFT(self.domain, norm=self._norm, workers=self._workers) class DistributedIFFT(DistributedFFTBase): @@ -149,9 +154,14 @@ def __init__(self, space, norm=None, workers=os.environ.get('OMP_NUM_THREADS', N assert isinstance(space, StencilVectorSpace) assert np.dtype(space.dtype).kind == 'c' workers = int(workers) if workers is not None else None + self._workers = workers + self._norm = norm super().__init__(space, lambda out: scifft.ifft( out, axis=1, overwrite_x=True, workers=workers, norm=norm)) + + def copy(self): + return DistributedIFFT(self.domain, norm=self._norm, workers=self._workers) class DistributedDCT(DistributedFFTBase): """ @@ -176,8 +186,14 @@ class DistributedDCT(DistributedFFTBase): """ def __init__(self, space, norm=None, workers=os.environ.get('OMP_NUM_THREADS', None), ttype=2): workers = int(workers) if workers is not None else None + self._workers = workers + self._norm = norm + self._ttype = ttype super().__init__(space, lambda out: scifft.dct( out, axis=1, overwrite_x=True, workers=workers, norm=norm, type=ttype)) + + def copy(self): + return DistributedDCT(self.domain, norm=self._norm, workers=self._workers, ttype=self._ttype) class DistributedIDCT(DistributedFFTBase): """ @@ -202,8 +218,14 @@ class DistributedIDCT(DistributedFFTBase): """ def __init__(self, space, norm=None, workers=os.environ.get('OMP_NUM_THREADS', None), ttype=2): workers = int(workers) if workers is not None else None + self._workers = workers + self._norm = norm + self._ttype = ttype super().__init__(space, lambda out: scifft.idct( out, axis=1, overwrite_x=True, workers=workers, norm=norm, type=ttype)) + + def copy(self): + return DistributedIDCT(self.domain, norm=self._norm, workers=self._workers, ttype=self._ttype) class DistributedDST(DistributedFFTBase): """ @@ -228,8 +250,14 @@ class DistributedDST(DistributedFFTBase): """ def __init__(self, space, norm=None, workers=os.environ.get('OMP_NUM_THREADS', None), ttype=2): workers = int(workers) if workers is not None else None + self._workers = workers + self._norm = norm + self._ttype = ttype super().__init__(space, lambda out: scifft.dst( out, axis=1, overwrite_x=True, workers=workers, norm=norm, type=ttype)) + + def copy(self): + return DistributedDST(self.domain, norm=self._norm, workers=self._workers, ttype=self._ttype) class DistributedIDST(DistributedFFTBase): """ @@ -254,5 +282,11 @@ class DistributedIDST(DistributedFFTBase): """ def __init__(self, space, norm=None, workers=os.environ.get('OMP_NUM_THREADS', None), ttype=2): workers = int(workers) if workers is not None else None + self._workers = workers + self._norm = norm + self._ttype = ttype super().__init__(space, lambda out: scifft.idst( out, axis=1, overwrite_x=True, workers=workers, norm=norm, type=ttype)) + + def copy(self): + return DistributedIDST(self.domain, norm=self._norm, workers=self._workers, ttype=self._ttype) diff --git a/psydac/linalg/kron.py b/psydac/linalg/kron.py index e4abfdd33..967dc963a 100644 --- a/psydac/linalg/kron.py +++ b/psydac/linalg/kron.py @@ -72,6 +72,9 @@ def ndim( self ): @property def mats( self ): return self._mats + + def copy(self): + return KroneckerStencilMatrix(self.domain, self.codomain, *self.mats) # ... def dot(self, x, out=None): @@ -257,6 +260,7 @@ def __init__(self, V, W, *args , with_pads=False): self._codomain = W self._mats = list(args) self._ndim = len(args) + self._with_pads = with_pads #-------------------------------------- # Abstract interface @@ -284,6 +288,9 @@ def ndim(self): @property def mats(self): return self._mats + + def copy(self): + return KroneckerDenseMatrix(self.domain, self.codomain, *self.mats, self._with_pads) # ... def dot(self, x, out=None): @@ -428,6 +435,9 @@ def __init__(self, V, W, solvers): # for now: allocate temporary arrays here (can be removed later) self._temp1, self._temp2 = self._allocate_temps() + + def copy(self): + return KroneckerLinearSolver(self.domain, self.codomain, self.solvers) def _setup_solvers(self): """ diff --git a/psydac/linalg/solvers.py b/psydac/linalg/solvers.py index 8c6505ca7..0f4d72bc3 100644 --- a/psydac/linalg/solvers.py +++ b/psydac/linalg/solvers.py @@ -131,6 +131,11 @@ def __init__(self, A, *, x0=None, tol=1e-6, maxiter=1000, verbose=False, recycle self._tmps = {key: self.domain.zeros() for key in ("v", "r", "p")} self._info = None + def copy(self): + return ConjugateGradient(self.linop, x0=self.get_options("x0"), tol=self.get_options("tol"), + maxiter=self.get_options("maxiter"), verbose=self.get_options("verbose"), + recycle=self.get_options("recycle")) + def solve(self, b, out=None): """ Conjugate gradient algorithm for solving linear system Ax=b. @@ -295,6 +300,11 @@ def __init__(self, A, *, pc=None, x0=None, tol=1e-6, maxiter=1000, verbose=False self._tmps = {**tmps_codomain, **tmps_domain} self._info = None + def copy(self): + return PConjugateGradient(self.linop, pc=self.get_options("pc"), x0=self.get_options("x0"), + tol=self.get_options("tol"), maxiter=self.get_options("maxiter"), + verbose=self.get_options("verbose"), recycle=self.get_options("recycle")) + def solve(self, b, out=None): """ Preconditioned Conjugate Gradient (PCG) solves the symetric positive definte @@ -454,6 +464,11 @@ def __init__(self, A, *, x0=None, tol=1e-6, maxiter=1000, verbose=False, recycle self._tmps = {key: self.domain.zeros() for key in ("v", "r", "p", "vs", "rs", "ps")} self._info = None + def copy(self): + return BiConjugateGradient(self.linop, x0=self.get_options("x0"), tol=self.get_options("tol"), + maxiter=self.get_options("maxiter"), verbose=self.get_options("verbose"), + recycle=self.get_options("recycle")) + def solve(self, b, out=None): """ Biconjugate gradient (BCG) algorithm for solving linear system Ax=b. @@ -643,6 +658,11 @@ def __init__(self, A, *, x0=None, tol=1e-6, maxiter=1000, verbose=False, recycle self._tmps = {key: self.domain.zeros() for key in ("v", "r", "p", "vr", "r0")} self._info = None + def copy(self): + return BiConjugateGradientStabilized(self.linop, x0=self.get_options("x0"), tol=self.get_options("tol"), + maxiter=self.get_options("maxiter"), verbose=self.get_options("verbose"), + recycle=self.get_options("recycle")) + def solve(self, b, out=None): """ Biconjugate gradient stabilized method (BCGSTAB) algorithm for solving linear system Ax=b. @@ -843,6 +863,11 @@ def __init__(self, A, *, pc=None, x0=None, tol=1e-6, maxiter=1000, verbose=False "rp0")} self._info = None + def copy(self): + return PBiConjugateGradientStabilized(self.linop, pc=self.get_options("pc"), x0=self.get_options("x0"), + tol=self.get_options("tol"), maxiter=self.get_options("maxiter"), + verbose=self.get_options("verbose"), recycle=self.get_options("recycle")) + def solve(self, b, out=None): """ Preconditioned biconjugate gradient stabilized method (PBCGSTAB) algorithm for solving linear system Ax=b. @@ -1085,6 +1110,11 @@ def __init__(self, A, *, x0=None, tol=1e-6, maxiter=1000, verbose=False, recycle self._tmps = {key: self.domain.zeros() for key in ("res_old", "res_new", "w_new", "w_work", "w_old", "v", "y")} self._info = None + def copy(self): + return MinimumResidual(self.linop, x0=self.get_options("x0"), tol=self.get_options("tol"), + maxiter=self.get_options("maxiter"), verbose=self.get_options("verbose"), + recycle=self.get_options("recycle")) + def solve(self, b, out=None): """ Use MINimum RESidual iteration to solve Ax=b @@ -1399,6 +1429,11 @@ def __init__(self, A, *, x0=None, tol=None, atol=None, btol=None, maxiter=1000, tmps_codomain = {key: self.codomain.zeros() for key in ("v", "v_work", "h", "hbar")} self._tmps = {**tmps_codomain, **tmps_domain} + def copy(self): + return LSMR(self.linop, x0=self.get_options("x0"), tol=self.get_options("tol"), atol=self.get_options("atol"), + btol=self.get_options("btol"), maxiter=self.get_options("maxiter"), conlim=self.get_options("conlim"), + verbose=self.get_options("verbose"), recycle=self.get_options("recycle")) + def get_success(self): return self._successful @@ -1740,6 +1775,11 @@ def __init__(self, A, *, x0=None, tol=1e-6, maxiter=100, verbose=False, recycle= self._Q = [] self._info = None + def copy(self): + return GMRES(self.linop, x0=self.get_options("x0"), tol=self.get_options("tol"), + maxiter=self.get_options("maxiter"), verbose=self.get_options("verbose"), + recycle=self.get_options("recycle")) + def solve(self, b, out=None): """ Generalized minimal residual algorithm for solving linear system Ax=b. @@ -1911,4 +1951,3 @@ def apply_givens_rotation(self, k, sn, cn): def dot(self, b, out=None): return self.solve(b, out=out) - diff --git a/psydac/linalg/tests/test_linalg.py b/psydac/linalg/tests/test_linalg.py index 56c826df1..dbd3fe70f 100644 --- a/psydac/linalg/tests/test_linalg.py +++ b/psydac/linalg/tests/test_linalg.py @@ -1,6 +1,12 @@ import pytest import numpy as np +from sympde.calculus import dot +from sympde.expr import BilinearForm, integral +from sympde.topology import Derham, elements_of, Square + +from psydac.api.discretization import discretize +from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL from psydac.linalg.block import BlockLinearOperator, BlockVector, BlockVectorSpace from psydac.linalg.basic import LinearOperator, ZeroOperator, IdentityOperator, ComposedLinearOperator, SumLinearOperator, PowerLinearOperator, ScaledLinearOperator from psydac.linalg.stencil import StencilVectorSpace, StencilVector, StencilMatrix @@ -52,7 +58,7 @@ def get_StencilVectorSpace(n1, n2, p1, p2, P1, P2): def get_positive_definite_StencilMatrix(V): - np.random.seed(2) + rng = np.random.default_rng(seed=2) assert isinstance(V, StencilVectorSpace) [n1, n2] = V._npts [p1, p2] = V._pads @@ -64,12 +70,12 @@ def get_positive_definite_StencilMatrix(V): for i in range(0, p1+1): if i != 0: for j in range(-p2, p2+1): - S[:, :, i, j] = 2*np.random.random()-1 + S[:, :, i, j] = 2*rng.random()-1 else: for j in range(1, p2+1): - S[:, :, i, j] = 2*np.random.random()-1 + S[:, :, i, j] = 2*rng.random()-1 S += S.T - S[:, :, 0, 0] = ((n1 * n2) - 1) / np.random.random() + S[:, :, 0, 0] = ((n1 * n2) - 1) / rng.random() S /= S[0, 0, 0, 0] S.remove_spurious_entries() @@ -627,7 +633,8 @@ def test_inverse_transpose_interaction(n1, n2, p1, p2, P1=False, P2=False): ### # Square root test - scaled_matrix = B * np.random.random() # Ensure the diagonal elements != 1 + rng = np.random.default_rng() + scaled_matrix = B * rng.random() # Ensure the diagonal elements != 1 diagonal_values = scaled_matrix.diagonal(sqrt=False).toarray() sqrt_diagonal_values = scaled_matrix.diagonal(sqrt=True).toarray() assert np.array_equal(sqrt_diagonal_values, np.sqrt(diagonal_values)) @@ -961,7 +968,7 @@ def test_internal_storage(): assert np.array_equal( y2_1.toarray(), y2_2.toarray() ) & np.array_equal( y2_2.toarray(), y2_3.toarray() ) #=============================================================================== -@pytest.mark.parametrize('solver', ['cg', 'pcg', 'bicg', 'minres', 'lsmr']) +@pytest.mark.parametrize('solver', ['cg', 'pcg', 'bicg', 'bicgstab', 'pbicgstab', 'minres', 'lsmr', 'gmres']) def test_x0update(solver): n1 = 4 @@ -980,7 +987,7 @@ def test_x0update(solver): # Create Inverse tol = 1e-6 - if solver == 'pcg': + if solver in ('pcg', 'pbicgstab'): A_inv = inverse(A, solver, pc=A.diagonal(inverse=True), tol=tol) else: A_inv = inverse(A, solver, tol=tol) @@ -1005,6 +1012,158 @@ def test_x0update(solver): x = A_inv.dot(b, out=b) assert A_inv.get_options('x0') is x +#=============================================================================== +@pytest.mark.parametrize('solver', ['cg', 'pcg', 'bicg', 'bicgstab', 'pbicgstab', 'minres', 'lsmr', 'gmres']) + +def test_ILO_copy(solver): + + rng = np.random.default_rng(seed=2) + domain = Square() + derham = Derham(domain, sequence=['h1', 'hcurl', 'l2']) + + domain_h = discretize(domain, ncells=[5,5], periodic=[False, False], comm=None) + derham_h = discretize(derham, domain_h, degree=[2, 2]) + + u, v = elements_of(derham.V1, names='u, v') + m1 = BilinearForm((u, v), integral(domain, dot(u, v))) + M1 = discretize(m1, domain_h, (derham_h.V1, derham_h.V1), backend=PSYDAC_BACKEND_GPYCCEL).assemble() + + if solver in ('pcg', 'pbicgstab'): + M1_inv = inverse(M1, solver, pc=M1.diagonal(inverse=True), tol=1e-3) + else: + M1_inv = inverse(M1, solver, tol=1e-3) + + M1_inv_copy = M1_inv.copy() + assert M1_inv is not M1_inv_copy + assert M1_inv is M1_inv # sanity check + + x = derham_h.V1.coeff_space.zeros() + rng.random(out=x[0]._data) + rng.random(out=x[1]._data) + + y1 = (M1_inv @ x).toarray() + y2 = (M1_inv_copy @ x).toarray() + assert np.allclose(y1, y2) + +#=============================================================================== +def test_LO_copy(): + + rng = np.random.default_rng(seed=2) + domain = Square() + derham = Derham(domain, sequence=['h1', 'hcurl', 'l2']) + + domain_h = discretize(domain, ncells=[5,5], periodic=[False, False], comm=None) + derham_h = discretize(derham, domain_h, degree=[2, 2]) + + u, v = elements_of(derham.V1, names='u, v') + m1 = BilinearForm((u, v), integral(domain, dot(u, v))) + M1 = discretize(m1, domain_h, (derham_h.V1, derham_h.V1), backend=PSYDAC_BACKEND_GPYCCEL).assemble() + I = IdentityOperator(derham_h.V1.coeff_space) + A = M1 + I + dt = 0.1 + + A1 = A # SumLinearOperator + A2 = A @ A # ComposedLinearOperator + A3 = A**3 # PowerLinearOperator + A4 = dt*A # ScaledLinearOperator + + B1 = A1.copy() + B2 = A2.copy() + B3 = A3.copy() + B4 = A4.copy() + + assert A1 is A1 # sanity + assert A1 is not B1 + assert A2 is not B2 + assert A3 is not B3 + assert A4 is not B4 + + x = derham_h.V1.coeff_space.zeros() + rng.random(out=x[0]._data) + rng.random(out=x[1]._data) + + assert np.allclose((A1@x).toarray(), (B1@x).toarray()) + assert np.allclose((A2@x).toarray(), (B2@x).toarray()) + assert np.allclose((A3@x).toarray(), (B3@x).toarray()) + assert np.allclose((A4@x).toarray(), (B4@x).toarray()) + +#=============================================================================== +def test_set_scalar(): + + n1 = 4 + n2 = 3 + p1 = 5 + p2 = 2 + P1 = False + P2 = False + V = get_StencilVectorSpace(n1, n2, p1, p2, P1, P2) + A = get_positive_definite_StencilMatrix(V) + I = IdentityOperator(V) + + x = V.zeros() + rng = np.random.default_rng() + rng.random(out=x._data) + + dt1 = 0.1 + dt2 = 0.2 + + A = A + dt1*I + y1 = A @ x + A.addends[1].set_scalar(dt2) + y2 = A @ x + diff = y2 - y1 + assert np.allclose(diff.toarray(), (0.1*x).toarray()) + +#=============================================================================== +def test_failing_BLO_add_cases(): + + domain = Square() + derham = Derham(domain, sequence=['h1', 'hcurl', 'l2']) + + domain_h = discretize(domain, ncells=[5,5], periodic=[False, False], comm=None) + derham_h = discretize(derham, domain_h, degree=[2, 2]) + + u, v = elements_of(derham.V1, names='u, v') + m1 = BilinearForm((u, v), integral(domain, dot(u, v))) + M1 = discretize(m1, domain_h, (derham_h.V1, derham_h.V1), backend=PSYDAC_BACKEND_GPYCCEL).assemble() + + V = derham_h.V1.coeff_space + + Z1 = ZeroOperator(V, V) + + A = M1 + A2 = BlockLinearOperator(V, V, ((None, None), (None, M1[1, 1]))) + + # 1. 'ZeroOperator' object has no attribute 'set_backend' + B = Z1 + S = A + B + + # 2. 'SumLinearOperator' object has no attribute 'set_backend' + B00 = IdentityOperator(V[0]) + B01 = StencilMatrix(V[1], V[0]) + B10 = StencilMatrix(V[0], V[1]) + + rng = np.random.default_rng() + rng.random(out=B01._data) + rng.random(out=B10._data) + + B[0, 0] = B00 + B[0, 1] = B01 + B[1, 0] = B10 + + S = A + B + + # 3. 'SumLinearOperator' object has no attribute 'copy' + B = BlockLinearOperator(V, V, ((M1[0,0] + IdentityOperator(V[0]), None), (None, None))) + + S = A2 + B + + # 4. 'ConjugateGradient' object has no attribute 'copy' + inv = inverse(M1[0, 0], 'cg', tol=1e-6, maxiter=100) + B = BlockLinearOperator(V, V, ((inv, None), (None, None))) + + S = A2 + B + #=============================================================================== # SCRIPT FUNCTIONALITY #=============================================================================== diff --git a/psydac/linalg/tests/test_linear_operator_rules.py b/psydac/linalg/tests/test_linear_operator_rules.py new file mode 100644 index 000000000..69de34217 --- /dev/null +++ b/psydac/linalg/tests/test_linear_operator_rules.py @@ -0,0 +1,98 @@ +from psydac.linalg.basic import IdentityOperator, SumLinearOperator, ScaledLinearOperator, ComposedLinearOperator +from psydac.linalg.block import BlockLinearOperator +from psydac.linalg.stencil import StencilMatrix + +from .test_linalg import get_StencilVectorSpace + + +def get_Hcurl_mass_matrix_2d(nc=5, comm=None): + + from sympde.calculus import dot + from sympde.expr import BilinearForm, integral + from sympde.topology import Derham, elements_of, Square + + from psydac.api.discretization import discretize + from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL + + domain = Square() + derham = Derham(domain, sequence=['h1', 'hcurl', 'l2']) + + domain_h = discretize(domain, ncells=[nc, nc], periodic=[False, False], comm=comm) + derham_h = discretize(derham, domain_h, degree=[2, 2]) + + u, v = elements_of(derham.V1, names='u, v') + m1 = BilinearForm((u, v), integral(domain, dot(u, v))) + M1 = discretize(m1, domain_h, (derham_h.V1, derham_h.V1), backend=PSYDAC_BACKEND_GPYCCEL).assemble() + + return M1 + + +def test_types_and_refs(): + + V = get_StencilVectorSpace(10, 5, p1=3, p2=2, P1=True, P2=False) + + I = IdentityOperator(V) # Immutable + S = StencilMatrix(V, V) # Mutable + P = StencilMatrix(V, V) # Mutable + + # TODO: set some entries of S, P to non-zero values + + # Example 1 + # --------- + + # Create simple sum + M = I + S + a = M + # New object, with references to [I, S] + # Type is SumLinearOperator + # If S is changed, so is M + assert isinstance(M, SumLinearOperator) + + # += does not modify the object, but creates a new one + M += 2*P # M := I + S + 2*P = a + 2*P + b = M + # New object, with references to the same atoms [I, S, P] + assert isinstance(M, SumLinearOperator) + assert M is not a + + # Store reference to M1 + M *= 2 # M := 2 * (I + S + 2*P) = 2 * b + # New object, with references to the same atoms [I, S, P] + assert isinstance(M, ScaledLinearOperator) + assert M is not b + + # Think about this one... are we OK with creating a new StencilMatrix? + W = S + 3 * S + P + assert isinstance(W, StencilMatrix) # currently... + + X = S @ S + assert isinstance(X, ComposedLinearOperator) + + assert isinstance(I @ S, StencilMatrix) + assert isinstance(S @ I, StencilMatrix) + assert (I @ S) is S + assert (S @ I) is S + + # Example 2 + # --------- + + M1 = get_Hcurl_mass_matrix_2d() + V1 = M1.domain + A = BlockLinearOperator(V1, V1, ((None, None), (None, M1[1, 1]))) + B = BlockLinearOperator(V1, V1, ((M1[0, 0] + IdentityOperator(V1[0]), None), (None, None))) + J = IdentityOperator(V1) + + # Sum: should we get a SumLinearOperator or a BlockLinearOperator? + C = A + B + r1 = C + assert isinstance(C, BlockLinearOperator) # currently... + + # In-place multiplication + C *= 5 # We want a new object here! + assert C is not r1 + assert isinstance(C, BlockLinearOperator) # debatable + + assert isinstance(J @ A, BlockLinearOperator) + assert isinstance(A @ J, BlockLinearOperator) + assert (J @ A) is A + assert (A @ J) is A diff --git a/psydac/linalg/tests/test_matrix_free.py b/psydac/linalg/tests/test_matrix_free.py index e475fb15f..1837b2a69 100644 --- a/psydac/linalg/tests/test_matrix_free.py +++ b/psydac/linalg/tests/test_matrix_free.py @@ -116,8 +116,14 @@ def test_solvers_matrix_free(solver): # Apply inverse and check y = A_inv @ x - error = np.linalg.norm( (b - y).toarray()) - assert np.linalg.norm( (b - y).toarray() ) < tol + x_approx = A @ y + error = np.linalg.norm((x - x_approx).toarray()) + if solver == 'lsmr': + # The LSMR solver has various termination conditions. In this case, + # convergence in the sence of error < tol is not achieved. + assert error < 10*tol + else: + assert error < tol #=============================================================================== # SCRIPT FUNCTIONALITY diff --git a/psydac/polar/dense.py b/psydac/polar/dense.py index 302e0a529..7c7cb0cdc 100644 --- a/psydac/polar/dense.py +++ b/psydac/polar/dense.py @@ -366,6 +366,9 @@ def domain(self): @property def codomain(self): return self._codomain + + def copy(self): + return DenseMatrix(self.domain, self.codomain, self._data.copy()) def transpose(self, conjugate=False): raise NotImplementedError()