Implement as_source_array/as_range_array for sparse NumpyMatrixOperators

src/pymor/operators/interface.py

@@ -8,6 +8,7 @@
 
 from pymor.algorithms import genericsolvers
 from pymor.core.base import abstractmethod
+from pymor.core.defaults import defaults
 from pymor.core.exceptions import InversionError, LinAlgError
 from pymor.parameters.base import ParametricObject
 from pymor.parameters.functionals import ParameterFunctional
@@ -368,6 +369,7 @@ def as_range_array(self, mu=None):
             The |VectorArray| defined above.
         """
         assert isinstance(self.source, NumpyVectorSpace) and self.linear
+        assert self.source.dim <= as_array_max_length()
         return self.apply(self.source.from_numpy(np.eye(self.source.dim)), mu=mu)
 
     def as_source_array(self, mu=None):
@@ -394,6 +396,7 @@ def as_source_array(self, mu=None):
             The |VectorArray| defined above.
         """
         assert isinstance(self.range, NumpyVectorSpace) and self.linear
+        assert self.range.dim <= as_array_max_length()
         return self.apply_adjoint(self.range.from_numpy(np.eye(self.range.dim)), mu=mu)
 
     def as_vector(self, mu=None):
@@ -541,7 +544,7 @@ def _add_sub(self, other, sign):
 
     def _radd_sub(self, other, sign):
         if other == 0:
-             return self
+            return self
         if not isinstance(other, Operator):
             return NotImplemented
 
@@ -584,3 +587,8 @@ def __neg__(self):
     def __str__(self):
         return f'{self.name}: R^{self.source.dim} --> R^{self.range.dim}  ' \
                f'(parameters: {self.parameters}, class: {self.__class__.__name__})'
+
+
+@defaults('value')
+def as_array_max_length(value=100):
+    return value

src/pymor/operators/numpy.py

@@ -221,12 +221,12 @@ def assemble(self, mu=None):
 
     def as_range_array(self, mu=None):
         if self.sparse:
-            raise NotImplementedError
+            return Operator.as_range_array(self)
         return self.range.from_numpy(self.matrix.T.copy())
 
     def as_source_array(self, mu=None):
         if self.sparse:
-            raise NotImplementedError
+            return Operator.as_source_array(self)
         return self.source.from_numpy(self.matrix.copy()).conj()
 
     def apply(self, U, mu=None):

src/pymortests/operators.py

@@ -13,6 +13,7 @@
 from pymor.operators.constructions import (SelectionOperator, InverseOperator, InverseAdjointOperator, IdentityOperator,
                                            LincombOperator, VectorArrayOperator)
 from pymor.operators.numpy import NumpyMatrixOperator
+from pymor.operators.interface import as_array_max_length
 from pymor.parameters.functionals import GenericParameterFunctional, ExpressionParameterFunctional
 from pymor.vectorarrays.block import BlockVectorSpace
 from pymor.vectorarrays.numpy import NumpyVectorSpace
@@ -484,3 +485,13 @@ def test_adjoint_vectorarray_op_apply_inverse_lstsq():
     v = V.to_numpy()
     u = np.linalg.lstsq(O, v.ravel())[0]
     assert np.all(almost_equal(U, U.space.from_numpy(u)))
+
+
+def test_as_range_array(operator_with_arrays):
+    op, mu, U, V = operator_with_arrays
+    if (not op.linear
+            or not isinstance(op.source, NumpyVectorSpace)
+            or op.source.dim > as_array_max_length()):
+        return
+    array = op.as_range_array(mu)
+    assert np.all(almost_equal(array.lincomb(U.to_numpy()), op.apply(U, mu=mu)))