From 31aa52b5c079a897b63fdbd9d1a5e77fd4b17136 Mon Sep 17 00:00:00 2001
From: jowezarek <julian.owezarek@tum.de>
Date: Thu, 13 Mar 2025 15:56:34 +0100
Subject: [PATCH 01/24] Add copy(), set_scalar and fix BLO __add__

---
 psydac/linalg/basic.py   | 26 ++++++++++++++++++++++++-
 psydac/linalg/block.py   |  3 ++-
 psydac/linalg/solvers.py | 41 +++++++++++++++++++++++++++++++++++++++-
 3 files changed, 67 insertions(+), 3 deletions(-)

diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py
index db28bfd7a..d455ff788 100644
--- a/psydac/linalg/basic.py
+++ b/psydac/linalg/basic.py
@@ -248,6 +248,10 @@ def codomain(self):
     def dtype(self):
         pass
 
+    @abstractmethod
+    def copy(self):
+        pass
+
     @abstractmethod
     def tosparse(self):
         pass
@@ -449,7 +453,7 @@ def __add__(self, B):
         assert isinstance(B, LinearOperator)
         assert self.domain == B.domain
         assert self.codomain == B.codomain
-        return B
+        return B.copy()
 
     def __sub__(self, B):
         assert isinstance(B, LinearOperator)
@@ -584,6 +588,14 @@ def operator(self):
     @property
     def dtype(self):
         return None
+    
+    def copy(self):
+        return ScaledLinearOperator(self.domain, self.codomain, self.scalar, self.operator)
+    
+    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() 
@@ -667,6 +679,9 @@ def addends(self):
     @property
     def dtype(self):
         return None
+    
+    def copy(self):
+        return SumLinearOperator(self.domain, self.codomain, *self.addends)
 
     def toarray(self):
         out = np.zeros(self.shape, dtype=self.dtype)
@@ -798,6 +813,9 @@ def multiplicants(self):
     @property
     def dtype(self):
         return None
+    
+    def copy(self):
+        return ComposedLinearOperator(self.domain, self.codomain, *self.multiplicants)
 
     def toarray(self):
         raise NotImplementedError('toarray() is not defined for ComposedLinearOperators.')
@@ -905,6 +923,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, self.factorial)
 
     def toarray(self):
         raise NotImplementedError('toarray() is not defined for PowerLinearOperators.')
@@ -1169,6 +1190,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 49e29b60d..3536b1446 100644
--- a/psydac/linalg/block.py
+++ b/psydac/linalg/block.py
@@ -1254,7 +1254,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/solvers.py b/psydac/linalg/solvers.py
index 3f81f64fa..5be05fee9 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)
-

From 2c0b9a79d68b67fab7c4664b9499623e9a760950 Mon Sep 17 00:00:00 2001
From: jowezarek <julian.owezarek@tum.de>
Date: Thu, 13 Mar 2025 16:51:05 +0100
Subject: [PATCH 02/24] add tests and copy to KroneckerLinearSolver

---
 psydac/linalg/basic.py             |   2 +-
 psydac/linalg/kron.py              |   3 +
 psydac/linalg/tests/test_linalg.py | 168 ++++++++++++++++++++++++++++-
 3 files changed, 168 insertions(+), 5 deletions(-)

diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py
index d455ff788..2949b3460 100644
--- a/psydac/linalg/basic.py
+++ b/psydac/linalg/basic.py
@@ -324,7 +324,7 @@ def __add__(self, B):
         """ Creates an object of the class SumLinearOperator unless B is a ZeroOperator in which case self is returned. """
         assert isinstance(B, LinearOperator)
         if isinstance(B, ZeroOperator):
-            return self
+            return self.copy()
         else:
             return SumLinearOperator(self.domain, self.codomain, self, B)
 
diff --git a/psydac/linalg/kron.py b/psydac/linalg/kron.py
index e4abfdd33..bd631bb13 100644
--- a/psydac/linalg/kron.py
+++ b/psydac/linalg/kron.py
@@ -428,6 +428,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/tests/test_linalg.py b/psydac/linalg/tests/test_linalg.py
index 3e03e8185..725e99a8b 100644
--- a/psydac/linalg/tests/test_linalg.py
+++ b/psydac/linalg/tests/test_linalg.py
@@ -7,6 +7,17 @@
 from psydac.linalg.solvers import ConjugateGradient, inverse
 from psydac.ddm.cart       import DomainDecomposition, CartDecomposition
 
+from sympde.calculus                import dot
+from sympde.expr                    import BilinearForm, integral
+from sympde.topology                import Cube, Derham, elements_of, Square
+
+from psydac.api.discretization      import discretize
+from psydac.api.settings            import PSYDAC_BACKENDS, PSYDAC_BACKEND_GPYCCEL
+from psydac.linalg.basic            import IdentityOperator, ZeroOperator
+from psydac.linalg.block            import BlockLinearOperator
+from psydac.linalg.solvers          import inverse
+from psydac.linalg.stencil          import StencilMatrix
+
 #===============================================================================
 
 n1array = [2, 7]
@@ -243,8 +254,9 @@ def test_square_stencil_basic(n1, n2, p1, p2, P1=False, P2=False):
     ## ___Addition and Substraction with ZeroOperators___
 
     # Adding a ZeroOperator does not change the StencilMatrix
-    assert (S + Z) is S
-    assert (Z + S) is S
+    # Update: We now return a copy now
+    #assert (S + Z) is S
+    #assert (Z + S) is S
 
     # Substracting a ZeroOperator and substracting from a ZeroOperator work as intended
     assert (S - Z) is S
@@ -961,7 +973,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 +992,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 +1017,154 @@ 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):
+
+    np.random.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.vector_space.zeros()
+    x[0]._data = np.random.random(x[0]._data.shape)
+    x[1]._data = np.random.random(x[1]._data.shape)
+
+    y1 = (M1_inv @ x).toarray()
+    y2 = (M1_inv_copy @ x).toarray()
+    assert np.allclose(y1, y2)
+
+#===============================================================================
+def test_LO_copy():
+
+    np.random.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.vector_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.vector_space.zeros()
+    x[0]._data = np.random.random(x[0]._data.shape)
+    x[1]._data = np.random.random(x[1]._data.shape)
+
+    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()
+    x._data = np.random.random(x._data.shape)
+    
+    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.vector_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])
+    B01._data = np.random.random(B01._data.shape)
+    B10 = StencilMatrix(V[0], V[1])
+    B10._data = np.random.random(B10._data.shape)
+    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
 #===============================================================================

From 9511c9dd99b1e567a2eab4d1d26c59d819163dd1 Mon Sep 17 00:00:00 2001
From: jowezarek <julian.owezarek@tum.de>
Date: Thu, 13 Mar 2025 18:39:12 +0100
Subject: [PATCH 03/24] add more copy() functions

---
 psydac/feec/multipatch/fem_linear_operators.py | 12 ++++++++++++
 psydac/feec/multipatch/operators.py            | 15 +++++++++++++++
 psydac/linalg/kron.py                          |  7 +++++++
 psydac/polar/dense.py                          |  3 +++
 4 files changed, 37 insertions(+)

diff --git a/psydac/feec/multipatch/fem_linear_operators.py b/psydac/feec/multipatch/fem_linear_operators.py
index 4b8da2553..1a57dec4b 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):
@@ -183,6 +189,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 +214,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 87368a1eb..dbdfb7af9 100644
--- a/psydac/feec/multipatch/operators.py
+++ b/psydac/feec/multipatch/operators.py
@@ -1011,6 +1011,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 +1101,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 +1122,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 +1141,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 +1160,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/kron.py b/psydac/linalg/kron.py
index bd631bb13..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):
diff --git a/psydac/polar/dense.py b/psydac/polar/dense.py
index 2e209f95e..e96b7ae1b 100644
--- a/psydac/polar/dense.py
+++ b/psydac/polar/dense.py
@@ -329,6 +329,9 @@ def domain(self):
     @property
     def codomain(self):
         return self._codomain
+    
+    def copy(self):
+        return DenseMatrix(self.domain, self.codomain, self._data)
 
     def transpose(self, conjugate=False):
         raise NotImplementedError()

From 4c2195b742611edd6d0d44fdf62da4c813b291c5 Mon Sep 17 00:00:00 2001
From: jowezarek <julian.owezarek@tum.de>
Date: Thu, 13 Mar 2025 19:57:08 +0100
Subject: [PATCH 04/24] add one more copy()

---
 psydac/linalg/fft.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/psydac/linalg/fft.py b/psydac/linalg/fft.py
index 5d67d178b..8edb87b67 100644
--- a/psydac/linalg/fft.py
+++ b/psydac/linalg/fft.py
@@ -26,6 +26,9 @@ def toarray(self):
 
     def tosparse(self):
         raise NotImplementedError('tosparse() is not defined for DistributedFFTBase.')
+    
+    def copy(self):
+        raise NotImplementedError('copy() is not implemented for DistributedFFTBase.')
 
     # Possible additions for the future:
     # * split off the LinearSolver class when used with the space ndarray (as used in the KroneckerLinearSolver),

From 52b71156d124f46c5b280ab5e6d5d272acccf447 Mon Sep 17 00:00:00 2001
From: jowezarek <julian.owezarek@tum.de>
Date: Fri, 14 Mar 2025 11:00:37 +0100
Subject: [PATCH 05/24] undo some ZeroOperator changes, two more instances of
 copy()

---
 psydac/feec/multipatch/fem_linear_operators.py |  3 +++
 psydac/feec/multipatch/operators.py            |  6 ++++++
 psydac/linalg/basic.py                         |  4 ++--
 psydac/linalg/block.py                         | 14 ++++++++------
 psydac/linalg/tests/test_linalg.py             |  5 ++---
 5 files changed, 21 insertions(+), 11 deletions(-)

diff --git a/psydac/feec/multipatch/fem_linear_operators.py b/psydac/feec/multipatch/fem_linear_operators.py
index 1a57dec4b..e9cbcbf0a 100644
--- a/psydac/feec/multipatch/fem_linear_operators.py
+++ b/psydac/feec/multipatch/fem_linear_operators.py
@@ -166,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 )
 
diff --git a/psydac/feec/multipatch/operators.py b/psydac/feec/multipatch/operators.py
index dbdfb7af9..935ce3034 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
diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py
index 2949b3460..4bbdca172 100644
--- a/psydac/linalg/basic.py
+++ b/psydac/linalg/basic.py
@@ -324,7 +324,7 @@ def __add__(self, B):
         """ Creates an object of the class SumLinearOperator unless B is a ZeroOperator in which case self is returned. """
         assert isinstance(B, LinearOperator)
         if isinstance(B, ZeroOperator):
-            return self.copy()
+            return self
         else:
             return SumLinearOperator(self.domain, self.codomain, self, B)
 
@@ -453,7 +453,7 @@ def __add__(self, B):
         assert isinstance(B, LinearOperator)
         assert self.domain == B.domain
         assert self.codomain == B.codomain
-        return B.copy()
+        return B
 
     def __sub__(self, B):
         assert isinstance(B, LinearOperator)
diff --git a/psydac/linalg/block.py b/psydac/linalg/block.py
index 3536b1446..ce1735572 100644
--- a/psydac/linalg/block.py
+++ b/psydac/linalg/block.py
@@ -728,6 +728,7 @@ def __mul__(self, a):
 
     # ...
     def __add__(self, M):
+        from psydac.linalg.basic import ZeroOperator
         if not isinstance(M, BlockLinearOperator):
             return LinearOperator.__add__(self, M)
 
@@ -737,9 +738,9 @@ 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()
-            else            : blocks[ij] = Bij + Mij
+            if   Bij is None or isinstance(Bij, ZeroOperator): blocks[ij] = Mij.copy() if Mij is not None else None
+            elif Mij is None or isinstance(Mij, ZeroOperator): blocks[ij] = Bij.copy()
+            else                                             : blocks[ij] = Bij + Mij
         mat = BlockLinearOperator(self.domain, self.codomain, blocks=blocks)
         if len(mat._blocks) != len(self._blocks):
             mat.set_backend(self._backend)
@@ -752,6 +753,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)
 
@@ -761,9 +763,9 @@ def __sub__(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
-            elif Mij is None: blocks[ij] =  Bij.copy()
-            else            : blocks[ij] =  Bij - Mij
+            if   Bij is None                                 : blocks[ij] = -Mij
+            elif Mij is None or isinstance(Mij, ZeroOperator): blocks[ij] =  Bij.copy()
+            else                                             : blocks[ij] =  Bij - Mij
         mat = BlockLinearOperator(self.domain, self.codomain, blocks=blocks)
         if len(mat._blocks) != len(self._blocks):
             mat.set_backend(self._backend)
diff --git a/psydac/linalg/tests/test_linalg.py b/psydac/linalg/tests/test_linalg.py
index 725e99a8b..ddebd2405 100644
--- a/psydac/linalg/tests/test_linalg.py
+++ b/psydac/linalg/tests/test_linalg.py
@@ -254,9 +254,8 @@ def test_square_stencil_basic(n1, n2, p1, p2, P1=False, P2=False):
     ## ___Addition and Substraction with ZeroOperators___
 
     # Adding a ZeroOperator does not change the StencilMatrix
-    # Update: We now return a copy now
-    #assert (S + Z) is S
-    #assert (Z + S) is S
+    assert (S + Z) is S
+    assert (Z + S) is S
 
     # Substracting a ZeroOperator and substracting from a ZeroOperator work as intended
     assert (S - Z) is S

From ec3d2c0dfcba05eef5baeebb976991d8b5442750 Mon Sep 17 00:00:00 2001
From: jowezarek <julian.owezarek@tum.de>
Date: Fri, 14 Mar 2025 11:02:38 +0100
Subject: [PATCH 06/24] fix imports in test_linalg

---
 psydac/linalg/tests/test_linalg.py | 17 ++++++-----------
 1 file changed, 6 insertions(+), 11 deletions(-)

diff --git a/psydac/linalg/tests/test_linalg.py b/psydac/linalg/tests/test_linalg.py
index ddebd2405..912e1e925 100644
--- a/psydac/linalg/tests/test_linalg.py
+++ b/psydac/linalg/tests/test_linalg.py
@@ -1,23 +1,18 @@
 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
 from psydac.linalg.solvers import ConjugateGradient, inverse
 from psydac.ddm.cart       import DomainDecomposition, CartDecomposition
 
-from sympde.calculus                import dot
-from sympde.expr                    import BilinearForm, integral
-from sympde.topology                import Cube, Derham, elements_of, Square
-
-from psydac.api.discretization      import discretize
-from psydac.api.settings            import PSYDAC_BACKENDS, PSYDAC_BACKEND_GPYCCEL
-from psydac.linalg.basic            import IdentityOperator, ZeroOperator
-from psydac.linalg.block            import BlockLinearOperator
-from psydac.linalg.solvers          import inverse
-from psydac.linalg.stencil          import StencilMatrix
-
 #===============================================================================
 
 n1array = [2, 7]

From 486aa3d8dc84568c3f7e8fb9407e9e76f7fe4971 Mon Sep 17 00:00:00 2001
From: jowezarek <julian.owezarek@tum.de>
Date: Thu, 20 Mar 2025 14:27:35 +0100
Subject: [PATCH 07/24] add copy() to 6 DistributedFFTBase subclasses

---
 psydac/linalg/fft.py | 37 ++++++++++++++++++++++++++++++++++---
 1 file changed, 34 insertions(+), 3 deletions(-)

diff --git a/psydac/linalg/fft.py b/psydac/linalg/fft.py
index 8edb87b67..dbe87a695 100644
--- a/psydac/linalg/fft.py
+++ b/psydac/linalg/fft.py
@@ -26,9 +26,6 @@ def toarray(self):
 
     def tosparse(self):
         raise NotImplementedError('tosparse() is not defined for DistributedFFTBase.')
-    
-    def copy(self):
-        raise NotImplementedError('copy() is not implemented for DistributedFFTBase.')
 
     # Possible additions for the future:
     # * split off the LinearSolver class when used with the space ndarray (as used in the KroneckerLinearSolver),
@@ -124,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):
@@ -152,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):
     """
@@ -179,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):
     """
@@ -205,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):
     """
@@ -231,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):
     """
@@ -257,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)

From 5191fcca70424803350f06bf5c895dd9d6a16831 Mon Sep 17 00:00:00 2001
From: jowezarek <julian.owezarek@tum.de>
Date: Thu, 20 Mar 2025 14:47:47 +0100
Subject: [PATCH 08/24] use np.random.default_rng and adjust copy method

---
 psydac/linalg/basic.py             | 12 +++++------
 psydac/linalg/tests/test_linalg.py | 33 +++++++++++++++++-------------
 psydac/polar/dense.py              |  2 +-
 3 files changed, 26 insertions(+), 21 deletions(-)

diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py
index 4bbdca172..a7c396558 100644
--- a/psydac/linalg/basic.py
+++ b/psydac/linalg/basic.py
@@ -275,8 +275,6 @@ def transpose(self, conjugate=False):
         """
         pass
 
-    # TODO: check if we should add a copy method!!!
-
     #-------------------------------------
     # Magic methods
     #-------------------------------------
@@ -590,7 +588,7 @@ def dtype(self):
         return None
     
     def copy(self):
-        return ScaledLinearOperator(self.domain, self.codomain, self.scalar, self.operator)
+        return ScaledLinearOperator(self.domain, self.codomain, self.scalar, self.operator.copy())
     
     def set_scalar(self, c):
         assert np.isscalar(c)
@@ -681,7 +679,8 @@ def dtype(self):
         return None
     
     def copy(self):
-        return SumLinearOperator(self.domain, self.codomain, *self.addends)
+        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)
@@ -815,7 +814,8 @@ def dtype(self):
         return None
     
     def copy(self):
-        return ComposedLinearOperator(self.domain, self.codomain, *self.multiplicants)
+        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.')
@@ -925,7 +925,7 @@ def factorial(self):
         return self._factorial
     
     def copy(self):
-        return PowerLinearOperator(self.domain, self.codomain, self.operator, self.factorial)
+        return PowerLinearOperator(self.domain, self.codomain, self.operator.copy(), self.factorial)
 
     def toarray(self):
         raise NotImplementedError('toarray() is not defined for PowerLinearOperators.')
diff --git a/psydac/linalg/tests/test_linalg.py b/psydac/linalg/tests/test_linalg.py
index 912e1e925..5e59f0307 100644
--- a/psydac/linalg/tests/test_linalg.py
+++ b/psydac/linalg/tests/test_linalg.py
@@ -58,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
@@ -70,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()
 
@@ -633,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))
@@ -1016,7 +1017,7 @@ def test_x0update(solver):
 
 def test_ILO_copy(solver):
 
-    np.random.seed(2)
+    rng = np.random.default_rng(seed=2)
     domain = Square()
     derham = Derham(domain, sequence=['h1', 'hcurl', 'l2'])
 
@@ -1037,8 +1038,8 @@ def test_ILO_copy(solver):
     assert M1_inv is M1_inv # sanity check
 
     x = derham_h.V1.vector_space.zeros()
-    x[0]._data = np.random.random(x[0]._data.shape)
-    x[1]._data = np.random.random(x[1]._data.shape)
+    rng.random(out=x[0]._data)
+    rng.random(out=x[1]._data)
 
     y1 = (M1_inv @ x).toarray()
     y2 = (M1_inv_copy @ x).toarray()
@@ -1047,7 +1048,7 @@ def test_ILO_copy(solver):
 #===============================================================================
 def test_LO_copy():
 
-    np.random.seed(2)
+    rng = np.random.default_rng(seed=2)
     domain = Square()
     derham = Derham(domain, sequence=['h1', 'hcurl', 'l2'])
 
@@ -1078,8 +1079,8 @@ def test_LO_copy():
     assert A4 is not B4
 
     x = derham_h.V1.vector_space.zeros()
-    x[0]._data = np.random.random(x[0]._data.shape)
-    x[1]._data = np.random.random(x[1]._data.shape)
+    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())
@@ -1100,7 +1101,8 @@ def test_set_scalar():
     I = IdentityOperator(V)
 
     x = V.zeros()
-    x._data = np.random.random(x._data.shape)
+    rng = np.random.default_rng()
+    rng.random(out=x._data)
     
     dt1 = 0.1
     dt2 = 0.2
@@ -1139,9 +1141,12 @@ def test_failing_BLO_add_cases():
     # 2. 'SumLinearOperator' object has no attribute 'set_backend'
     B00 = IdentityOperator(V[0])
     B01 = StencilMatrix(V[1], V[0])
-    B01._data = np.random.random(B01._data.shape)
     B10 = StencilMatrix(V[0], V[1])
-    B10._data = np.random.random(B10._data.shape)
+
+    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
diff --git a/psydac/polar/dense.py b/psydac/polar/dense.py
index e96b7ae1b..6f9edf82c 100644
--- a/psydac/polar/dense.py
+++ b/psydac/polar/dense.py
@@ -331,7 +331,7 @@ def codomain(self):
         return self._codomain
     
     def copy(self):
-        return DenseMatrix(self.domain, self.codomain, self._data)
+        return DenseMatrix(self.domain, self.codomain, self._data.copy())
 
     def transpose(self, conjugate=False):
         raise NotImplementedError()

From a6d758b5cc818468c8a70dfef83541c6a7f3453d Mon Sep 17 00:00:00 2001
From: jowezarek <julian.owezarek@tum.de>
Date: Thu, 20 Mar 2025 14:58:20 +0100
Subject: [PATCH 09/24] remove three copy() instances in BLO

---
 psydac/linalg/block.py | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/psydac/linalg/block.py b/psydac/linalg/block.py
index ce1735572..4b4c32311 100644
--- a/psydac/linalg/block.py
+++ b/psydac/linalg/block.py
@@ -738,9 +738,9 @@ 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 or isinstance(Bij, ZeroOperator): blocks[ij] = Mij.copy() if Mij is not None else None
-            elif Mij is None or isinstance(Mij, ZeroOperator): blocks[ij] = Bij.copy()
-            else                                             : blocks[ij] = Bij + Mij
+            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):
             mat.set_backend(self._backend)
@@ -763,9 +763,9 @@ def __sub__(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
-            elif Mij is None or isinstance(Mij, ZeroOperator): blocks[ij] =  Bij.copy()
-            else                                             : blocks[ij] =  Bij - Mij
+            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):
             mat.set_backend(self._backend)
@@ -994,7 +994,7 @@ def __iadd__(self, M):
 
             Bij = self[ij]
             if Bij is None:
-                self[ij] = Mij.copy()
+                self[ij] = Mij
             else:
                 Bij += Mij
 

From 644455df99e135f2db5b02994e0488b423f7f5e3 Mon Sep 17 00:00:00 2001
From: jowezarek <julian.owezarek@tum.de>
Date: Fri, 21 Mar 2025 18:01:05 +0100
Subject: [PATCH 10/24] re-add copy() to __iadd__ method

---
 psydac/linalg/block.py | 9 ++++++++-
 1 file changed, 8 insertions(+), 1 deletion(-)

diff --git a/psydac/linalg/block.py b/psydac/linalg/block.py
index 4b4c32311..5357f5af8 100644
--- a/psydac/linalg/block.py
+++ b/psydac/linalg/block.py
@@ -994,7 +994,14 @@ def __iadd__(self, M):
 
             Bij = self[ij]
             if Bij is None:
-                self[ij] = Mij
+                # 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
 

From fb60defdbbfd48d1dfe971083efbd2c03c95c0d5 Mon Sep 17 00:00:00 2001
From: jowezarek <julian.owezarek@tum.de>
Date: Tue, 25 Mar 2025 11:34:58 +0100
Subject: [PATCH 11/24] remove ZeroLO import

---
 psydac/linalg/block.py | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/psydac/linalg/block.py b/psydac/linalg/block.py
index 5357f5af8..1bfa34395 100644
--- a/psydac/linalg/block.py
+++ b/psydac/linalg/block.py
@@ -728,7 +728,6 @@ def __mul__(self, a):
 
     # ...
     def __add__(self, M):
-        from psydac.linalg.basic import ZeroOperator
         if not isinstance(M, BlockLinearOperator):
             return LinearOperator.__add__(self, M)
 
@@ -738,8 +737,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
-            elif Mij is None: blocks[ij] = Bij
+            if   Bij is None: blocks[ij] = Mij.copy()
+            elif Mij is None: blocks[ij] = Bij.copy()
             else            : blocks[ij] = Bij + Mij
         mat = BlockLinearOperator(self.domain, self.codomain, blocks=blocks)
         if len(mat._blocks) != len(self._blocks):
@@ -764,7 +763,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
+            elif Mij is None: blocks[ij] =  Bij.copy()
             else            : blocks[ij] =  Bij - Mij
         mat = BlockLinearOperator(self.domain, self.codomain, blocks=blocks)
         if len(mat._blocks) != len(self._blocks):

From 580f70d6d84438daa2e15e6d79f9464aefd752c2 Mon Sep 17 00:00:00 2001
From: jowezarek <julian.owezarek@tum.de>
Date: Tue, 25 Mar 2025 13:49:43 +0100
Subject: [PATCH 12/24] fix faulty test

---
 psydac/linalg/tests/test_matrix_free.py | 10 ++++++++--
 1 file changed, 8 insertions(+), 2 deletions(-)

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

From 55264fce5ae701d634d523bd1bf47e3d908a89b9 Mon Sep 17 00:00:00 2001
From: jowezarek <julian.owezarek@tum.de>
Date: Tue, 25 Mar 2025 13:50:50 +0100
Subject: [PATCH 13/24] re-undo copy() instances in block __add__ and __sub__

---
 psydac/linalg/block.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/psydac/linalg/block.py b/psydac/linalg/block.py
index 1bfa34395..112083729 100644
--- a/psydac/linalg/block.py
+++ b/psydac/linalg/block.py
@@ -737,8 +737,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):
@@ -763,7 +763,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):

From 44d6475f018d94435bc571f64dad3a27c64bde95 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= <yaman.guclu@gmail.com>
Date: Wed, 16 Apr 2025 15:17:40 +0200
Subject: [PATCH 14/24] Create test_linear_operator_rules.py

First attempt at defining the rules for operating with `LinearOperator` objects.
---
 linalg/tests/test_linear_operator_rules.py | 48 ++++++++++++++++++++++
 1 file changed, 48 insertions(+)
 create mode 100644 linalg/tests/test_linear_operator_rules.py

diff --git a/linalg/tests/test_linear_operator_rules.py b/linalg/tests/test_linear_operator_rules.py
new file mode 100644
index 000000000..2561d395d
--- /dev/null
+++ b/linalg/tests/test_linear_operator_rules.py
@@ -0,0 +1,48 @@
+from psydac.linalg.basic   import IdentityOperator, SumLinearOperator, ScaledLinearOperator, ComposedLinearOperator
+#from psydac.linalg.block   import BlockVectorSpace
+from psydac.linalg.stencil import StencilVectorSpace, StencilMatrix
+
+from .test_linalg import get_StencilVectorSpace
+
+
+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
+  M1 = I + S
+  a = M1
+  # New object, with references to [I, S]
+  # Type is SumLinearOperator
+  # If S is changed, so is M1
+  assert isinstance(M1, SumLinearOperator)
+
+  # 
+  M1 += 2*P  # M1 := I + S + 2*P
+  b = M1
+  # Same object, with references to [I, S, P]
+  assert isinstance(M1, SumLinearOperator)
+  assert a is b
+
+  # Store reference to M1
+  M1 *= 2 # -> 2 * (I + S + 2*P)
+  # New object, with references to [I, S, P]
+  # Not the same object as b := I + S + 2*P
+  assert isinstance(M1, ScaledLinearOperator)
+  assert b is not M1
+
+  # Think about this one...
+  M2 = S + S
+  assert isinstance(M2, StencilMatrix)
+
+  M3 = S @ S
+  assert isinstance(M3, ComposedLinearOperator)

From 586b6b7d387f38dc65a06a199eff410d3f866914 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= <yaman.guclu@gmail.com>
Date: Wed, 16 Apr 2025 20:29:14 +0200
Subject: [PATCH 15/24] Update test_linalg.py

Replace `vector_space` with `coeff_space`.
---
 psydac/linalg/tests/test_linalg.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/psydac/linalg/tests/test_linalg.py b/psydac/linalg/tests/test_linalg.py
index 5e59f0307..21a80e382 100644
--- a/psydac/linalg/tests/test_linalg.py
+++ b/psydac/linalg/tests/test_linalg.py
@@ -1037,7 +1037,7 @@ def test_ILO_copy(solver):
     assert M1_inv is not M1_inv_copy
     assert M1_inv is M1_inv # sanity check
 
-    x = derham_h.V1.vector_space.zeros()
+    x = derham_h.V1.coeff_space.zeros()
     rng.random(out=x[0]._data)
     rng.random(out=x[1]._data)
 
@@ -1058,7 +1058,7 @@ def test_LO_copy():
     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.vector_space)
+    I  = IdentityOperator(derham_h.V1.coeff_space)
     A = M1 + I
     dt = 0.1
 
@@ -1078,7 +1078,7 @@ def test_LO_copy():
     assert A3 is not B3
     assert A4 is not B4
 
-    x = derham_h.V1.vector_space.zeros()
+    x = derham_h.V1.coeff_space.zeros()
     rng.random(out=x[0]._data)
     rng.random(out=x[1]._data)
 
@@ -1127,7 +1127,7 @@ def test_failing_BLO_add_cases():
     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.vector_space
+    V = derham_h.V1.coeff_space
 
     Z1 = ZeroOperator(V, V)
 

From ead08e1f6ab8b6e449ab960548ab64ba27e6b564 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= <yaman.guclu@gmail.com>
Date: Wed, 23 Apr 2025 14:13:06 +0200
Subject: [PATCH 16/24] Only run linear algebra tests

---
 .github/workflows/continuous-integration.yml | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/.github/workflows/continuous-integration.yml b/.github/workflows/continuous-integration.yml
index ae4c2e09d..379c00ec5 100644
--- a/.github/workflows/continuous-integration.yml
+++ b/.github/workflows/continuous-integration.yml
@@ -185,7 +185,7 @@ jobs:
         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()

From 6868e45f6d8732d430d86043f30d3d6b72508d81 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= <yaman.guclu@gmail.com>
Date: Wed, 23 Apr 2025 14:14:31 +0200
Subject: [PATCH 17/24] Move test_linear_operator_rules.py to the correct
 folder

---
 {linalg => psydac/linalg}/tests/test_linear_operator_rules.py | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 rename {linalg => psydac/linalg}/tests/test_linear_operator_rules.py (100%)

diff --git a/linalg/tests/test_linear_operator_rules.py b/psydac/linalg/tests/test_linear_operator_rules.py
similarity index 100%
rename from linalg/tests/test_linear_operator_rules.py
rename to psydac/linalg/tests/test_linear_operator_rules.py

From 43121cdb8cd0d07fa72acb761f8a9ad879b6bfd3 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= <yaman.guclu@gmail.com>
Date: Wed, 23 Apr 2025 14:59:00 +0200
Subject: [PATCH 18/24] Fix test_linear_operator_rules.py

---
 .../tests/test_linear_operator_rules.py       | 27 +++++++++++--------
 1 file changed, 16 insertions(+), 11 deletions(-)

diff --git a/psydac/linalg/tests/test_linear_operator_rules.py b/psydac/linalg/tests/test_linear_operator_rules.py
index 2561d395d..5e46e07b4 100644
--- a/psydac/linalg/tests/test_linear_operator_rules.py
+++ b/psydac/linalg/tests/test_linear_operator_rules.py
@@ -26,23 +26,28 @@ def test_types_and_refs():
   # If S is changed, so is M1
   assert isinstance(M1, SumLinearOperator)
 
-  # 
-  M1 += 2*P  # M1 := I + S + 2*P
+  # += does not modify the object, but creates a new one
+  M1 += 2*P  # M1 := I + S + 2*P = a + 2*P
   b = M1
-  # Same object, with references to [I, S, P]
+  # New object, with references to the same atoms [I, S, P]
   assert isinstance(M1, SumLinearOperator)
-  assert a is b
+  assert M1 is not a
 
   # Store reference to M1
-  M1 *= 2 # -> 2 * (I + S + 2*P)
-  # New object, with references to [I, S, P]
-  # Not the same object as b := I + S + 2*P
+  M1 *= 2   # M1 := 2 * (I + S + 2*P) = 2 * b
+  # New object, with references to the same atoms [I, S, P]
   assert isinstance(M1, ScaledLinearOperator)
-  assert b is not M1
+  assert M1 is not b
 
-  # Think about this one...
-  M2 = S + S
-  assert isinstance(M2, StencilMatrix)
+  # Think about this one... are we OK with creating a new StencilMatrix?
+  M2 = S + 3 * S + P
+  assert isinstance(M2, StencilMatrix)  # ?
 
   M3 = S @ S
   assert isinstance(M3, ComposedLinearOperator)
+
+  # Example 2
+  # ---------
+
+  
+  

From c78be77df2c8b27583ae7ae1fa332e28f93a3734 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= <yaman.guclu@gmail.com>
Date: Wed, 23 Apr 2025 15:21:10 +0200
Subject: [PATCH 19/24] Add BlockLinearOperator test in
 test_linear_operator_rules.py

---
 .../tests/test_linear_operator_rules.py       | 62 ++++++++++++++-----
 1 file changed, 45 insertions(+), 17 deletions(-)

diff --git a/psydac/linalg/tests/test_linear_operator_rules.py b/psydac/linalg/tests/test_linear_operator_rules.py
index 5e46e07b4..84f764ce2 100644
--- a/psydac/linalg/tests/test_linear_operator_rules.py
+++ b/psydac/linalg/tests/test_linear_operator_rules.py
@@ -1,10 +1,26 @@
 from psydac.linalg.basic   import IdentityOperator, SumLinearOperator, ScaledLinearOperator, ComposedLinearOperator
 #from psydac.linalg.block   import BlockVectorSpace
 from psydac.linalg.stencil import StencilVectorSpace, StencilMatrix
+from psydac.api.settings   import PSYDAC_BACKEND_GPYCCEL
 
 from .test_linalg import get_StencilVectorSpace
 
 
+def get_Hcurl_mass_matrix_2d(nc=5, comm=None):
+
+    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)
@@ -19,35 +35,47 @@ def test_types_and_refs():
   # ---------
 
   # Create simple sum
-  M1 = I + S
-  a = M1
+  M = I + S
+  a = M
   # New object, with references to [I, S]
   # Type is SumLinearOperator
-  # If S is changed, so is M1
-  assert isinstance(M1, SumLinearOperator)
+  # If S is changed, so is M
+  assert isinstance(M, SumLinearOperator)
 
   # += does not modify the object, but creates a new one
-  M1 += 2*P  # M1 := I + S + 2*P = a + 2*P
-  b = M1
+  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(M1, SumLinearOperator)
-  assert M1 is not a
+  assert isinstance(M, SumLinearOperator)
+  assert M is not a
 
   # Store reference to M1
-  M1 *= 2   # M1 := 2 * (I + S + 2*P) = 2 * b
+  M *= 2   # M := 2 * (I + S + 2*P) = 2 * b
   # New object, with references to the same atoms [I, S, P]
-  assert isinstance(M1, ScaledLinearOperator)
-  assert M1 is not b
+  assert isinstance(M, ScaledLinearOperator)
+  assert M is not b
 
   # Think about this one... are we OK with creating a new StencilMatrix?
-  M2 = S + 3 * S + P
-  assert isinstance(M2, StencilMatrix)  # ?
+  W = S + 3 * S + P
+  assert isinstance(W, StencilMatrix)  # currently...
 
-  M3 = S @ S
-  assert isinstance(M3, ComposedLinearOperator)
+  X = S @ S
+  assert isinstance(X, ComposedLinearOperator)
 
   # 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)))
+
+  # 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

From aa4846594a9c3ad2832024b8a6c30391e9fc5e3c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= <yaman.guclu@gmail.com>
Date: Wed, 23 Apr 2025 15:37:42 +0200
Subject: [PATCH 20/24] Remove __imul__, __iadd__, __isub__ from
 BlockLinearOperator

---
 psydac/linalg/block.py | 68 +++++++++++++++++++++---------------------
 1 file changed, 34 insertions(+), 34 deletions(-)

diff --git a/psydac/linalg/block.py b/psydac/linalg/block.py
index 112083729..9a2f23cb8 100644
--- a/psydac/linalg/block.py
+++ b/psydac/linalg/block.py
@@ -972,27 +972,27 @@ 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)
+    # 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
+    #     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:
+    #         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)
@@ -1000,33 +1000,33 @@ def __iadd__(self, 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
+    #             self[ij] = Mij.copy()
+    #         else:
+    #             Bij += Mij
 
-        return self
+    #     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):

From 4bc54bb73682fdaafa4bcf2b193d74ceadccc4d4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= <yaman.guclu@gmail.com>
Date: Wed, 23 Apr 2025 16:13:26 +0200
Subject: [PATCH 21/24] Add missing imports into test_linear_operator_rules.py

---
 psydac/linalg/tests/test_linear_operator_rules.py | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/psydac/linalg/tests/test_linear_operator_rules.py b/psydac/linalg/tests/test_linear_operator_rules.py
index 84f764ce2..e52dd3f7a 100644
--- a/psydac/linalg/tests/test_linear_operator_rules.py
+++ b/psydac/linalg/tests/test_linear_operator_rules.py
@@ -1,13 +1,19 @@
 from psydac.linalg.basic   import IdentityOperator, SumLinearOperator, ScaledLinearOperator, ComposedLinearOperator
 #from psydac.linalg.block   import BlockVectorSpace
 from psydac.linalg.stencil import StencilVectorSpace, StencilMatrix
-from psydac.api.settings   import PSYDAC_BACKEND_GPYCCEL
 
 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'])
 

From 4a2ab8a55c2fc8420ea17c63d0240a7f64ec4df9 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= <yaman.guclu@gmail.com>
Date: Wed, 23 Apr 2025 16:28:04 +0200
Subject: [PATCH 22/24] Add missing import in test_linear_operator_rules.py

---
 psydac/linalg/tests/test_linear_operator_rules.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/psydac/linalg/tests/test_linear_operator_rules.py b/psydac/linalg/tests/test_linear_operator_rules.py
index e52dd3f7a..e25846c3b 100644
--- a/psydac/linalg/tests/test_linear_operator_rules.py
+++ b/psydac/linalg/tests/test_linear_operator_rules.py
@@ -1,6 +1,6 @@
 from psydac.linalg.basic   import IdentityOperator, SumLinearOperator, ScaledLinearOperator, ComposedLinearOperator
-#from psydac.linalg.block   import BlockVectorSpace
-from psydac.linalg.stencil import StencilVectorSpace, StencilMatrix
+from psydac.linalg.block   import BlockLinearOperator
+from psydac.linalg.stencil import StencilMatrix
 
 from .test_linalg import get_StencilVectorSpace
 

From 68ca9067343b69342460b6c41d413a53225de4c3 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= <yaman.guclu@gmail.com>
Date: Wed, 23 Apr 2025 17:04:33 +0200
Subject: [PATCH 23/24] Run fewer tests on macOS

---
 .github/workflows/continuous-integration.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/continuous-integration.yml b/.github/workflows/continuous-integration.yml
index 379c00ec5..f9f7b9636 100644
--- a/.github/workflows/continuous-integration.yml
+++ b/.github/workflows/continuous-integration.yml
@@ -179,7 +179,7 @@ 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'

From c1fd8c8992a634862cfba67703b5999fde1260b0 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= <yaman.guclu@gmail.com>
Date: Wed, 7 May 2025 15:34:32 +0200
Subject: [PATCH 24/24] Add rules on composition with identity

---
 psydac/linalg/tests/test_linear_operator_rules.py | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/psydac/linalg/tests/test_linear_operator_rules.py b/psydac/linalg/tests/test_linear_operator_rules.py
index e25846c3b..69de34217 100644
--- a/psydac/linalg/tests/test_linear_operator_rules.py
+++ b/psydac/linalg/tests/test_linear_operator_rules.py
@@ -68,6 +68,11 @@ def test_types_and_refs():
   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
   # ---------
 
@@ -75,6 +80,7 @@ def test_types_and_refs():
   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
@@ -85,3 +91,8 @@ def test_types_and_refs():
   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