Repository: pymor/pymor

Issue:
BlockOperator.apply_inverse infinite recursion with a special LincombOperator as a block
Here is an example:

```python
from pymor.operators.block import BlockOperator
from pymor.operators.constructions import IdentityOperator
from pymor.vectorarrays.numpy import NumpyVectorSpace

I = IdentityOperator(NumpyVectorSpace(1))
B = BlockOperator([[I * (-1)]])
v = B.source.ones()

B.apply_inverse(v)
```

The issue is that `assemble` for `I * (-1)` always returns a new operator object.

Discussion:

Patch:
diff --git a/src/pymor/algorithms/lincomb.py b/src/pymor/algorithms/lincomb.py
index 8f1928878..fa827a736 100644
--- a/src/pymor/algorithms/lincomb.py
+++ b/src/pymor/algorithms/lincomb.py
@@ -26,11 +26,7 @@ def assemble_lincomb(operators, coefficients, solver_options=None, name=None):
     where `O_i` are |Operators| and `c_i` scalar coefficients.
 
     This function is called in the :meth:`assemble` method of |LincombOperator| and
-    is not intended to be used directly. The assembled |Operator| is expected to
-    no longer be a |LincombOperator| nor should it contain any |LincombOperators|.
-    If an assembly of the given linear combination is not possible, `None` is returned.
-    The special case of a |LincombOperator| with a single operator (i.e. a scaled |Operator|)
-    is allowed as assemble_lincomb implements `apply_inverse` for this special case.
+    is not intended to be used directly.
 
     To form the linear combination of backend |Operators| (containing actual matrix data),
     :meth:`~pymor.operators.interface.Operator._assemble_lincomb` will be called
@@ -49,7 +45,7 @@ def assemble_lincomb(operators, coefficients, solver_options=None, name=None):
 
     Returns
     -------
-    The assembled |Operator| if assembly is possible, otherwise `None`.
+    The assembled |Operator|.
     """
     return AssembleLincombRules(tuple(coefficients), solver_options, name).apply(tuple(operators))
 
@@ -289,10 +285,6 @@ class AssembleLincombRules(RuleTable):
             raise RuleNotMatchingError
         return op
 
-    @match_generic(lambda ops: len(ops) == 1)
-    def action_only_one_operator(self, ops):
-        return LincombOperator(ops, self.coefficients, name=self.name)
-
     @match_always
-    def action_failed(self, ops):
-        return None
+    def action_return_lincomb(self, ops):
+        return LincombOperator(ops, self.coefficients, name=self.name)
diff --git a/src/pymor/operators/constructions.py b/src/pymor/operators/constructions.py
index d7b22d161..f689128d8 100644
--- a/src/pymor/operators/constructions.py
+++ b/src/pymor/operators/constructions.py
@@ -138,31 +138,32 @@ class LincombOperator(Operator):
         from pymor.algorithms.lincomb import assemble_lincomb
         operators = tuple(op.assemble(mu) for op in self.operators)
         coefficients = self.evaluate_coefficients(mu)
+        # try to form a linear combination
         op = assemble_lincomb(operators, coefficients, solver_options=self.solver_options,
                               name=self.name + '_assembled')
-        if op:
+        # To avoid infinite recursions, only use the result if at least one of the following
+        # is true:
+        #   - The operator is parametric, so the the result of assemble *must* be a different,
+        #     non-parametric operator.
+        #   - One of self.operators changed by calling 'assemble' on it.
+        #   - The result of assemble_lincomb is of a different type than the original operator.
+        #   - assemble_lincomb could simplify the list of assembled operators,
+        #     which we define to be the case when the number of operators has ben reduced.
+        if (self.parametric
+                or operators != self.operators  # for this comparison to work self.operators always has to be a tuple!
+                or type(op) != type(self)
+                or len(op.operators) < len(operators)):
             return op
         else:
-            if self.parametric or operators != self.operators:
-                return LincombOperator(operators, coefficients, solver_options=self.solver_options,
-                                       name=self.name + '_assembled')
-            else:  # this can only happen when both operators and self.operators are tuples!
-                return self  # avoid infinite recursion
+            return self
 
     def jacobian(self, U, mu=None):
-        from pymor.algorithms.lincomb import assemble_lincomb
         if self.linear:
             return self.assemble(mu)
         jacobians = [op.jacobian(U, mu) for op in self.operators]
-        coefficients = self.evaluate_coefficients(mu)
         options = self.solver_options.get('jacobian') if self.solver_options else None
-        jac = assemble_lincomb(jacobians, coefficients, solver_options=options,
-                               name=self.name + '_jacobian')
-        if jac is None:
-            return LincombOperator(jacobians, coefficients, solver_options=options,
-                                   name=self.name + '_jacobian')
-        else:
-            return jac
+        return LincombOperator(jacobians, self.coefficients, solver_options=options,
+                               name=self.name + '_jacobian').assemble(mu)
 
     def d_mu(self, parameter, index=0):
         for op in self.operators:


Test patch:
diff --git a/src/pymortests/operators.py b/src/pymortests/operators.py
index e0a6f45d9..5c613a2c1 100644
--- a/src/pymortests/operators.py
+++ b/src/pymortests/operators.py
@@ -466,3 +466,15 @@ def test_as_range_array(operator_with_arrays):
         return
     array = op.as_range_array(mu)
     assert np.all(almost_equal(array.lincomb(U.to_numpy()), op.apply(U, mu=mu)))
+
+
+def test_issue_1276():
+    from pymor.operators.block import BlockOperator
+    from pymor.operators.constructions import IdentityOperator
+    from pymor.vectorarrays.numpy import NumpyVectorSpace
+
+    I = IdentityOperator(NumpyVectorSpace(1))
+    B = BlockOperator([[I * (-1)]])
+    v = B.source.ones()
+
+    B.apply_inverse(v)