diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py
index dc86957b4..33a0200e2 100644
--- a/psydac/linalg/basic.py
+++ b/psydac/linalg/basic.py
@@ -432,26 +432,13 @@ def __matmul__(self, B):
 #===============================================================================
 class IdentityOperator(LinearOperator):
 
-    def __new__(cls, domain, codomain=None):
+    def __init__(self, domain, codomain=None):
 
         assert isinstance(domain, VectorSpace)
         if codomain:
             assert isinstance(codomain, VectorSpace)
             assert domain == codomain
 
-        from psydac.linalg.block import BlockVectorSpace, BlockLinearOperator
-        if isinstance(domain, BlockVectorSpace):
-            spaces = domain.spaces
-            blocks = {}
-            for i, V in enumerate(spaces):
-                blocks[i,i] = IdentityOperator(V)
-            return BlockLinearOperator(domain, domain, blocks)
-        else:
-            return super().__new__(cls)
-
-    
-    def __init__(self, domain, codomain=None):
-
         self._domain = domain
         self._codomain = domain
 
diff --git a/psydac/linalg/tests/test_linalg.py b/psydac/linalg/tests/test_linalg.py
index d929e42cf..ea00cbf3e 100644
--- a/psydac/linalg/tests/test_linalg.py
+++ b/psydac/linalg/tests/test_linalg.py
@@ -368,12 +368,11 @@ def test_square_block_basic(n1, n2, p1, p2, P1=False, P2=False):
     assert isinstance(B - B1, BlockLinearOperator)
 
     # Adding and Substracting BlockLOs and other LOs returns a SumLinearOperator object
-    # Update 21.12.: ZeroLOs and IdentityLOs from and/or to BlockVectorSpaces are now BlockLOs
-    # thus the sums/differences below should be BlockLOs again
-    assert isinstance(B + BI, BlockLinearOperator)
-    assert isinstance(BI + B, BlockLinearOperator)
-    assert isinstance(B - BI, BlockLinearOperator)
-    assert isinstance(BI - B, BlockLinearOperator)
+    # reverted changes made in #261
+    assert isinstance(B + BI, SumLinearOperator)
+    assert isinstance(BI + B, SumLinearOperator)
+    assert isinstance(B - BI, SumLinearOperator)
+    assert isinstance(BI - B, SumLinearOperator)
 
     # Negating a BlockLO works as intended
     assert isinstance(-B, BlockLinearOperator)
@@ -426,19 +425,19 @@ def test_square_block_basic(n1, n2, p1, p2, P1=False, P2=False):
     assert isinstance(BZ @ B, ComposedLinearOperator)
 
     # Composing a BlockLO with the IdentityOperator does not change the object
-    # due to the 21.12. change not valid anymore
-    #assert B @ BI == B
-    #assert BI @ B == B
+    # reverted changes from #261
+    assert B @ BI == B
+    assert BI @ B == B
     # but: 
-    assert array_equal((B @ BI) @ vb, B @ vb)
-    assert array_equal((BI @ B) @ vb, B @ vb)
+    # assert array_equal((B @ BI) @ vb, B @ vb)
+    # assert array_equal((BI @ B) @ vb, B @ vb)
 
     ## ___Raising to the power of 0 and 1___
 
     # Raising a BlockLO to the power of 1 or 0 does not change the object / returns an IdentityOperator
-    # 21.12. change: B**0 a BlockLO with IdentityLOs at the diagonal
+    # reverted changes from #261
     assert B**1 is B
-    assert isinstance(B**0, BlockLinearOperator)
+    assert isinstance(B**0, IdentityOperator)
     assert sparse_equal(B**0, BI)
 
 #===============================================================================