Improve behavior of LinearOperator objects

[jowezarek]

Closes #453

The issue described in #453 is not a ZeroOperator bug, but rather an issue in the __add__ method of the BlockLinearOperator:

1     def __add__(self, M):
2         if not isinstance(M, BlockLinearOperator):
3             return LinearOperator.__add__(self, M)
4
5         assert M.  domain is self.domain
6         assert M.codomain is self.codomain
7         blocks  = {}
8         for ij in set(self._blocks.keys()) | set(M._blocks.keys()):
9             Bij = self[ij]
10            Mij = M[ij]
11            if   Bij is None: blocks[ij] = Mij.copy()
12            elif Mij is None: blocks[ij] = Bij.copy()
13            else            : blocks[ij] = Bij + Mij
14        mat = BlockLinearOperator(self.domain, self.codomain, blocks=blocks)
15        if len(mat._blocks) != len(self._blocks):
16            mat.set_backend(self._backend)
17        ...

First Issue

In lines 11 & 12, it is assumed that all LinearOperators have a copy method. This is a useful requirement, as we don't want the last line of the following code

A = BlockLinearOperator(V, V, ( (A1, None), (None, A2) )
B = BlockLinearOperator(V, V, ( (B1, None), (None, None) )
C = A + B                    # line 12 for ij = [1,1]: C[i, j] = A2.copy()
A2 *= 2

to change C.
In fact, this code does not run if A2 does not have an implemented copy method.
Hence I made copy a mandatory method for all LinearOperators.

Second Issue

Line 16 is called whenever there is a block [i,j] such that in the case C = A + B A[i,j] == None and B[i,j] != None.
The method set_backend looks like this:

1     def set_backend(self, backend):
2         if isinstance(self.domain, BlockVectorSpace) and isinstance(self.domain.spaces[0], BlockVectorSpace):
3             return
4
5         if isinstance(self.codomain, BlockVectorSpace) and isinstance(self.codomain.spaces[0], BlockVectorSpace):
6             return
7
8         if backend is None:return
9         if backend is self._backend:return
10
11        from psydac.api.ast.linalg import LinearOperatorDot
12        from psydac.linalg.stencil import StencilInterfaceMatrix, StencilMatrix
13
14        if not all(isinstance(b, (StencilMatrix, StencilInterfaceMatrix)) for b in self._blocks.values()):
15            for b in self._blocks.values():
16                b.set_backend(backend)
17            return
18            ...

Here, unless one of the four break criteria is met, and unless all not-None blocks are StencilMatrices or StencilInterfaceMatrices, each LinearOperator must have a set_backend function defined. I don't think that is useful, and hence replaced line 16 as follows

                  if getattr(b, "set_backend", None) is not None:
                      b.set_backend(backend)

Now, only if a set_backend method exists, it is called.

Third Issue

In line 13 of the BlockLinearOperator.__add__ method, if either one of the LinearOperators Bij or Mij is a ZeroOperator we have Bij + Mij = Mij or Bij + Mij = Bij rather than Mij.copy() or Bij.copy(). The same issue appears in the __sub__ method. Hence I changed the code as follows:

# in __add__
            if   Bij is None or isinstance(Bij, ZeroOperator): blocks[ij] = Mij.copy() if Mij is not None else None
            elif Mij is None or isinstance(Mij, ZeroOperator): blocks[ij] = Bij.copy()
            else                                             : blocks[ij] = Bij + Mij
# in __sub__
            if   Bij is None                                 : blocks[ij] = -Mij
            elif Mij is None or isinstance(Mij, ZeroOperator): blocks[ij] =  Bij.copy()
            else                                             : blocks[ij] =  Bij - Mij

Additional Change

I have required, and used locally, a set_scalar method in the ScaledLinearOperator class for a long time already.
It looks as follows and is included in this PR

    def set_scalar(self, c):
        assert np.isscalar(c)
        assert np.iscomplexobj(c) == (self.codomain._dtype == complex)
        self._scalar = c

Finally, I have added many instances of copy methods - sometimes like so

# in ScaledLinearOperator
    def copy(self):
        return ScaledLinearOperator(self.domain, self.codomain, self.scalar, self.operator)

and sometimes like so

# in DistributedFFTBase
    def copy(self):
        raise NotImplementedError('copy() is not implemented for DistributedFFTBase.')

Tests

are included in this PR for "standard" LinearOperator copy functions. I've also included tests that would previously fail, either due to a non-existent set_backend method or due to a non-existent copy method.

[codacy-production]

Coverage summary from Codacy

See diff coverage on Codacy

Coverage variation	Diff coverage
✅ -24.35%	✅ 64.58%

Coverage variation details

	Coverable lines	Covered lines	Coverage
Common ancestor commit (f96e888)	30771	18605	60.46%
Head commit (04208b7)	61662 (+30891)	22268 (+3663)	36.11% (-24.35%)

Coverage variation is the difference between the coverage for the head and common ancestor commits of the pull request branch: <coverage of head commit> - <coverage of common ancestor commit>

Diff coverage details

	Coverable lines	Covered lines	Diff coverage
Pull request (#473)	96	62	64.58%

Diff coverage is the percentage of lines that are covered by tests out of the coverable lines that the pull request added or modified: <covered lines added or modified>/<coverable lines added or modified> * 100%

See your quality gate settings    Change summary preferences

[yguclu]

Good job @jowezarek! Thanks!

I have a few comments and questions

[yguclu]

@jowezarek Is this PR ready for being reviewed one last time? If yes, please update the PR description with a list of the changes made, and the issues which it fixes. When merging with devel, I will use that part of the PR description to create the commit message.