Sparse BlockOperator

src/pymor/algorithms/lincomb.py

@@ -177,49 +177,29 @@ def scaled_iden_coeff(op):
         so_op = SecondOrderModelOperator(alpha, beta, A, B)
         return so_op
 
-    @match_class_all(BlockDiagonalOperator)
-    def action_BlockDiagonalOperator(self, ops):
-        coefficients = self.coefficients
-        num_source_blocks = ops[0].num_source_blocks
-        blocks = np.empty((num_source_blocks,), dtype=object)
-        if len(ops) > 1:
-            for i in range(num_source_blocks):
-                operators_i = [op.blocks[i, i] for op in ops]
-                blocks[i] = assemble_lincomb(operators_i, coefficients,
-                                             solver_options=self.solver_options, name=self.name)
-                if blocks[i] is None:
-                    return None
-            return BlockDiagonalOperator(blocks)
-        else:
-            c = coefficients[0]
-            if c == 1:
-                return ops[0]
-            for i in range(num_source_blocks):
-                blocks[i] = ops[0].blocks[i, i] * c
-            return BlockDiagonalOperator(blocks)
-
     @match_class_all(BlockOperatorBase)
     def action_BlockOperatorBase(self, ops):
         coefficients = self.coefficients
         shape = ops[0].blocks.shape
-        blocks = np.empty(shape, dtype=object)
+        blocks = np.zeros(shape, dtype=object)
         operator_type = ((BlockOperator if ops[0].blocked_source else BlockColumnOperator) if ops[0].blocked_range
                          else BlockRowOperator)
         if len(ops) > 1:
-            for (i, j) in np.ndindex(shape):
-                operators_ij = [op.blocks[i, j] for op in ops]
+            # the sparsity pattern can differ.
+            merged_coords = np.unique(np.vstack([op.block_coords for op in ops]), axis=0)
+            for (i, j) in merged_coords:
+                # currently we use a self-written slicing for coo_matrix in scipy.sparse
+                # TODO: find a better way to do it
+                operators_ij = [op.slice((i, j)) for op in ops if (i, j) in op.block_coords]
                 blocks[i, j] = assemble_lincomb(operators_ij, coefficients,
                                                 solver_options=self.solver_options, name=self.name)
-                if blocks[i, j] is None:
-                    return None
-            return operator_type(blocks)
         else:
             c = coefficients[0]
             if c == 1:
                 return ops[0]
-            for (i, j) in np.ndindex(shape):
-                blocks[i, j] = ops[0].blocks[i, j] * c
-            return operator_type(blocks)
+            for i, j, op in zip(ops[0].blocks.row, ops[0].blocks.col, ops[0].blocks.data):
+                blocks[i, j] = op * c
+        return operator_type(blocks)
 
     @match_generic(lambda ops: sum(1 for op in ops if isinstance(op, LowRankOperator)) >= 2)
     def action_merge_low_rank_operators(self, ops):

src/pymor/algorithms/projection.py

@@ -232,19 +232,28 @@ def action_BlockOperatorBase(self, op):
                 range_bases = self.range_basis.blocks
             else:
                 range_bases = [None] * len(op.range.subspaces)
+            range_dims = [len(rb) if rb is not None else op.range.subspaces[i].dim
+                          for i, rb in enumerate(range_bases)]
         else:
             range_bases = [self.range_basis]
+            range_dims = [op.range.dim]
         if op.blocked_source:
             if self.source_basis is not None:
                 source_bases = self.source_basis.blocks
             else:
                 source_bases = [None] * len(op.source.subspaces)
+            source_dims = [len(sb) if sb is not None else op.source.subspaces[i].dim
+                           for i, sb in enumerate(source_bases)]
         else:
             source_bases = [self.source_basis]
-
-        projected_ops = np.array([[project(op.blocks[i, j], rb, sb)
-                                   for j, sb in enumerate(source_bases)]
-                                  for i, rb in enumerate(range_bases)])
+            source_dims = [op.source.dim]
+
+        # NOTE: this does not use a potential sparsity pattern of the BlockOperator
+        # need source and range dims for the case where None is involved.
+        projected_ops = np.array([[project(op.slice((i, j)), rb, sb) if (i, j) in op.block_coords
+                                   else NumpyMatrixOperator(np.zeros((rdim, sdim)))
+                                   for j, (sb, sdim) in enumerate(zip(source_bases, source_dims))]
+                                  for i, (rb, rdim) in enumerate(zip(range_bases, range_dims))])
         if self.range_basis is None and op.blocked_range:
             return BlockColumnOperator(np.sum(projected_ops, axis=1))
         elif self.source_basis is None and op.blocked_source:

src/pymor/algorithms/to_matrix.py

@@ -95,15 +95,15 @@ def action_NumpyMatrixOperator(self, op):
     @match_class(BlockOperatorBase)
     def action_BlockOperator(self, op):
         format = self.format
-        op_blocks = op.blocks
+        # Note: to_dense() only puts zero operators for all None entries, is this needed?
+        op = op.to_dense()
         mat_blocks = [[] for i in range(op.num_range_blocks)]
         is_dense = True
-        for i in range(op.num_range_blocks):
-            for j in range(op.num_source_blocks):
-                mat_ij = self.apply(op_blocks[i, j])
-                if sps.issparse(mat_ij):
-                    is_dense = False
-                mat_blocks[i].append(mat_ij)
+        for i, j, op in zip(op.blocks.row, op.blocks.col, op.blocks.data):
+            mat_ij = self.apply(op)
+            if sps.issparse(mat_ij):
+                is_dense = False
+            mat_blocks[i].append(mat_ij)
         if format is None and is_dense or format == 'dense':
             return np.block(mat_blocks)
         else: