fix functionality for applying to a multidimensional domain

src/operators/matrix_product_operator.py

@@ -104,8 +104,8 @@ class MatrixProductOperator(EndomorphicOperator):
             appl_space_shape = tuple(appl_space_shape)
             self._active_axes = tuple(active_axes)
 
-            self._mat_last_n = tuple([-appl_dim + i for i in range(appl_dim)])
-            self._mat_first_n = np.arange(appl_dim)
+        self._mat_last_n = tuple([-appl_dim + i for i in range(appl_dim)])
+        self._mat_first_n = np.arange(appl_dim)
 
         # Test if the matrix and the array it will be applied to fit
         if matrix.shape[:appl_dim] != appl_space_shape:
@@ -126,7 +126,14 @@ class MatrixProductOperator(EndomorphicOperator):
 
         if self._spaces is None:
             if not self._flatten:
-                res = m.dot(x.val)
+                if times:
+                    res = np.tensordot(m,x.val,axes=len(x.shape))
+                else:
+                    mat_axes = np.flip(self._mat_last_n)
+                    field_axes = self._mat_first_n
+                    res = np.tensordot(m,x.val,axes=(mat_axes,field_axes))
+                    res = res.reshape(np.flip(res.shape))
+
             else:
                 res = m.dot(x.val.flatten()).reshape(self._domain.shape)
             return Field(self._domain, res)