Make KroneckerLinearSolver a subclass of LinearOperator

docs/source/modules/linalg/block.rst

@@ -4,7 +4,6 @@ linalg.block
    * :ref:`BlockVectorSpace <blockvectorspace>`
    * :ref:`BlockVector <blockvector>`
    * :ref:`BlockLinearOperator <blocklinearoperator>`
-   * :ref:`BlockDiagonalSolver <blockdiagonalsolver>`  
 
 .. inheritance-diagram:: psydac.linalg.block
 
@@ -31,11 +30,3 @@ BlockLinearOperator
 
 .. autoclass:: psydac.linalg.block.BlockLinearOperator
    :members:
-
-.. _blockdiagonalsolver:
-
-BlockDiagonalSolver
--------------------
-
-.. autoclass:: psydac.linalg.block.BlockDiagonalSolver
-   :members:

docs/source/modules/linalg/direct_solvers.rst

@@ -1,20 +1,11 @@
 linalg.direct_solvers
 =====================
 
-   * :ref:`DirectSolver <directsolver>`
    * :ref:`BandedSolver <bandedsolver>`
    * :ref:`SparseSolver <sparsesolver>`
 
 .. inheritance-diagram:: psydac.linalg.direct_solvers
 
-.. _directsolver:
-
-DirectSolver
-------------
-
-.. autoclass:: psydac.linalg.direct_solvers.DirectSolver
-   :members:
-
 .. _bandedsolver:
 
 BandedSolver

psydac/feec/global_projectors.py

@@ -4,7 +4,7 @@
 
 from psydac.linalg.kron           import KroneckerLinearSolver
 from psydac.linalg.stencil        import StencilVector
-from psydac.linalg.block          import BlockDiagonalSolver, BlockVector
+from psydac.linalg.block          import BlockLinearOperator, BlockVector
 from psydac.core.bsplines         import quadrature_grid
 from psydac.utilities.quadratures import gauss_legendre
 from psydac.fem.basic             import FemField
@@ -168,7 +168,7 @@ def __init__(self, space, nquads = None):
             self._grid_x += [block_x]
             self._grid_w += [block_w]
 
-            solverblocks += [KroneckerLinearSolver(tensorspaces[i].vector_space, solvercells)]
+            solverblocks += [KroneckerLinearSolver(tensorspaces[i].vector_space, tensorspaces[i].vector_space, solvercells)]
 
             dataslice = tuple(slice(p, -p) for p in tensorspaces[i].vector_space.pads)
             dofs[i] = rhsblocks[i]._data[dataslice]
@@ -177,11 +177,13 @@ def __init__(self, space, nquads = None):
         args = (*intp_x, *quad_x, *quad_w, *dofs)
         self._func = lambda *fun: func(*args, *fun)
 
-        # build a BlockDiagonalSolver, if necessary
+        # build a BlockLinearOperator, if necessary
         if len(solverblocks) == 1:
             self._solver = solverblocks[0]
         else:
-            self._solver = BlockDiagonalSolver(self._space.vector_space, blocks=solverblocks)
+            domain = codomain = self._space.vector_space
+            blocks = {(i, i): B_i for i, B_i in enumerate(solverblocks)}
+            self._solver = BlockLinearOperator(domain, codomain, blocks)
 
     @property
     def space(self):
@@ -307,7 +309,7 @@ def __call__(self, fun):
         else:
             self._func(fun)
 
-        coeffs = self._solver.solve(self._rhs)
+        coeffs = self._solver.dot(self._rhs)
 
         return FemField(self._space, coeffs=coeffs)
 

psydac/linalg/basic.py

@@ -1005,5 +1005,14 @@ def space(self):
         pass
 
     @abstractmethod
-    def solve(self, rhs, out=None, transposed=False):
+    def transpose(self):
+        """Return the transpose of the LinearSolver."""
         pass
+
+    @abstractmethod
+    def solve(self, rhs, out=None):
+        pass
+
+    @property
+    def T(self):
+        return self.transpose()

psydac/linalg/block.py

@@ -7,12 +7,12 @@
 from types import MappingProxyType
 from scipy.sparse import bmat, lil_matrix
 
-from psydac.linalg.basic  import VectorSpace, Vector, LinearOperator, LinearSolver, ZeroOperator
-from psydac.ddm.cart      import InterfaceCartDecomposition
-from psydac.ddm.utilities import get_data_exchanger
-from psydac.linalg.stencil import StencilVector, StencilMatrix
+from psydac.linalg.basic    import VectorSpace, Vector, LinearOperator
+from psydac.linalg.stencil  import StencilMatrix
+from psydac.ddm.cart        import InterfaceCartDecomposition
+from psydac.ddm.utilities   import get_data_exchanger
 
-__all__ = ('BlockVectorSpace', 'BlockVector', 'BlockLinearOperator', 'BlockDiagonalSolver')
+__all__ = ('BlockVectorSpace', 'BlockVector', 'BlockLinearOperator')
 
 #===============================================================================
 class BlockVectorSpace(VectorSpace):
@@ -1321,138 +1321,3 @@ def func(blocks, v, out, **args):
 
         self._func    = func
         self._backend = backend
-
-#===============================================================================
-class BlockDiagonalSolver( LinearSolver ):
-    """
-    A LinearSolver that can be written as blocks of other LinearSolvers,
-    i.e. it can be seen as a solver for linear equations with block-diagonal matrices.
-
-    The space of this solver has to be of the type BlockVectorSpace.
-
-    Parameters
-    ----------
-    V : psydac.linalg.block.BlockVectorSpace
-        Space of the new blocked linear solver.
-
-    blocks : dict | list | tuple
-        LinearSolver objects (optional).
-
-        a) 'blocks' can be dictionary with
-            . key   = integer i >= 0
-            . value = corresponding LinearSolver Lii
-
-        b) 'blocks' can be list of LinearSolvers (or tuple of these) where blocks[i]
-            is the LinearSolver Lii (if None, we assume null operator)
-
-    """
-    def __init__( self, V, blocks=None ):
-
-        assert isinstance( V, BlockVectorSpace )
-
-        self._space   = V
-        self._nblocks = V.n_blocks
-
-        # Store blocks in list (hence, they can be manually changed later)
-        self._blocks   = [None] * self._nblocks
-
-        if blocks:
-
-            if isinstance( blocks, dict ):
-                for i, L in blocks.items():
-                    self[i] = L
-
-            elif isinstance( blocks, (list, tuple) ):
-                for i, L in enumerate( blocks ):
-                    self[i] = L
-
-            else:
-                raise ValueError( "Blocks can only be given as dict or 1D list/tuple." )
-
-    #--------------------------------------
-    # Abstract interface
-    #--------------------------------------
-    @property
-    def space( self ):
-        """
-        The space this BlockDiagonalSolver works on.
-        """
-        return self._space
-
-    # ...
-    def solve( self, rhs, out=None, transposed=False ):
-        """
-        Solves the linear system for the given right-hand side rhs.
-        An out vector can be supplied, otherwise a new vector will be allocated.
-
-        This operation supports in-place operations, given that the underlying solvers
-        do as well.
-
-        Parameters
-        ----------
-        rhs : BlockVector
-            The input right-hand side.
-
-        out : BlockVector | NoneType
-            The output vector, or None.
-
-        transposed : Bool
-            If true, and supported by the underlying solvers,
-            rhs is solved against the transposed right-hand sides.
-
-        Returns
-        -------
-        out : BlockVector
-            Either `out`, if given as input parameter, or a newly-allocated vector.
-            In all cases, it holds the result of the computation.
-        """
-        assert isinstance(rhs, BlockVector)
-        assert rhs.space is self.space
-        if out is None:
-            out = self.space.zeros()
-
-        assert isinstance(out, BlockVector)
-        assert out.space is self.space
-
-        rhs.update_ghost_regions()
-
-        for i, L in enumerate(self._blocks):
-            if L is None:
-                raise NotImplementedError('All solvers have to be defined.')
-            L.solve(rhs[i], out=out[i], transposed=transposed)
-
-        return out
-
-    #--------------------------------------
-    # Other properties/methods
-    #--------------------------------------
-    @property
-    def blocks( self ):
-        """ Immutable 1D view (tuple) of the linear solvers,
-            including the empty blocks as 'None' objects.
-        """
-        return tuple(self._blocks)
-    # ...
-    @property
-    def n_blocks( self ):
-        """
-        The number of blocks in the matrix.
-        """
-        return self._nblocks
-
-    # ...
-    def __getitem__( self, key ):
-        assert 0 <= key < self._nblocks
-
-        return self._blocks.get( key, None )
-
-    # ...
-    def __setitem__( self, key, value ):
-        assert 0 <= key < self._nblocks
-
-        assert isinstance( value, LinearSolver )
-
-        # Check domain of rhs
-        assert value.space is self.space[key]
-
-        self._blocks[key] = value