diff --git a/psydac/api/tests/test_2d_multipatch_mapping_maxwell.py b/psydac/api/tests/test_2d_multipatch_mapping_maxwell.py
index 90dbc6b0a..a9759ba86 100644
--- a/psydac/api/tests/test_2d_multipatch_mapping_maxwell.py
+++ b/psydac/api/tests/test_2d_multipatch_mapping_maxwell.py
@@ -93,10 +93,19 @@ def run_maxwell_2d(uex, f, alpha, domain, *, ncells=None, degree=None, filename=
     equation_h = discretize(equation, domain_h, [Vh, Vh], backend=PSYDAC_BACKEND_GPYCCEL)
     l2norm_h   = discretize(l2norm, domain_h, Vh, backend=PSYDAC_BACKEND_GPYCCEL)
 
-    equation_h.set_solver('pcg', pc='jacobi', tol=1e-8)
+    # Explicitly assemble the linear system
+    equation_h.assemble()
 
+    # Extract the matrix' diagonal to build a Jacobi preconditioner
+    jacobi_pc = equation_h.linear_system.lhs.diagonal(inverse=True)
+
+    # Choose a linear solver and pass any flags to it
+    equation_h.set_solver('pcg', pc=jacobi_pc, tol=1e-8)
+
+    # Solve the linear system and obtain the solution as a FEM field
     uh = equation_h.solve()
 
+    # Compute the L2 norm of the error
     l2_error = l2norm_h.assemble(F=uh)
 
     return l2_error, uh
diff --git a/psydac/api/tests/test_api_2d_compatible_spaces.py b/psydac/api/tests/test_api_2d_compatible_spaces.py
index 1638af024..e84e3782e 100644
--- a/psydac/api/tests/test_api_2d_compatible_spaces.py
+++ b/psydac/api/tests/test_api_2d_compatible_spaces.py
@@ -399,7 +399,8 @@ def run_maxwell_time_harmonic_2d_dir(uex, f, alpha, ncells, degree):
     #+++++++++++++++++++++++++++++++
 
     # Solve linear system
-    M_inv = inverse(M, 'pcg', pc='jacobi', tol=1e-8)
+    jacobi_pc = M.diagonal(inverse=True)
+    M_inv = inverse(M, 'pcg', pc=jacobi_pc, tol=1e-8)
     sol = M_inv @ b
 
     uh       = FemField( Vh, sol )
diff --git a/psydac/api/tests/test_api_feec_1d.py b/psydac/api/tests/test_api_feec_1d.py
index 260bb0381..3c06204c7 100644
--- a/psydac/api/tests/test_api_feec_1d.py
+++ b/psydac/api/tests/test_api_feec_1d.py
@@ -336,7 +336,8 @@ def discrete_energies(self, e, b):
         step_ampere_1d = dt * ( M0_dir_inv @ D0_T_dir @ M1 )
 
     elif bc_mode == 'penalization':
-        M0_M0_bc_inv = inverse(M0+M0_bc, 'pcg', pc='jacobi', **kwargs)
+        M0_M0_bc = M0 + M0_bc
+        M0_M0_bc_inv = inverse(M0_M0_bc, 'pcg', pc = M0_M0_bc.diagonal(inverse=True), **kwargs)
         step_ampere_1d = dt * ( M0_M0_bc_inv @ D0_T @ M1 )
 
     half_step_faraday_1d = (dt/2) * D0
diff --git a/psydac/api/tests/test_api_feec_2d.py b/psydac/api/tests/test_api_feec_2d.py
index 5c27b678e..1453211f1 100644
--- a/psydac/api/tests/test_api_feec_2d.py
+++ b/psydac/api/tests/test_api_feec_2d.py
@@ -519,7 +519,8 @@ def discrete_energies(self, e, b):
         M1_inv = inverse(M1, 'cg', **kwargs)
         step_ampere_2d = dt * (M1_inv @ D1_T @ M2)
     else:
-        M1_M1_bc_inv = inverse(M1 + M1_bc, 'pcg', pc='jacobi', **kwargs)
+        M1_M1_bc = M1 + M1_bc
+        M1_M1_bc_inv = inverse(M1_M1_bc, 'pcg', pc = M1_M1_bc.diagonal(inverse=True), **kwargs)
         step_ampere_2d = dt * (M1_M1_bc_inv @ D1_T @ M2)
 
     half_step_faraday_2d = (dt/2) * D1
diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py
index 352421d0b..352b4e380 100644
--- a/psydac/linalg/basic.py
+++ b/psydac/linalg/basic.py
@@ -3,12 +3,26 @@
 # Copyright 2018 Yaman Güçlü, Jalal Lakhlili
 # Copyright 2022 Yaman Güçlü, Said Hadjout, Julian Owezarek
 
-from abc   import ABC, abstractmethod
-from scipy.sparse import coo_matrix
+from abc import ABC, abstractmethod
+
 import numpy as np
+from scipy.sparse import coo_matrix
 
-__all__ = ('VectorSpace', 'Vector', 'LinearOperator', 'ZeroOperator', 'IdentityOperator', 'ScaledLinearOperator',
-           'SumLinearOperator', 'ComposedLinearOperator', 'PowerLinearOperator', 'InverseLinearOperator', 'LinearSolver')
+from psydac.utilities.utils import is_real
+
+__all__ = (
+    'VectorSpace',
+    'Vector',
+    'LinearOperator',
+    'ZeroOperator',
+    'IdentityOperator',
+    'ScaledLinearOperator',
+    'SumLinearOperator',
+    'ComposedLinearOperator',
+    'PowerLinearOperator',
+    'InverseLinearOperator',
+    'LinearSolver'
+)
 
 #===============================================================================
 class VectorSpace(ABC):
@@ -448,7 +462,6 @@ def __new__(cls, domain, codomain=None):
             return BlockLinearOperator(domain, domain, blocks)
         else:
             return super().__new__(cls)
-
     
     def __init__(self, domain, codomain=None):
 
@@ -876,14 +889,51 @@ def dot(self, v, out=None):
 #===============================================================================
 class InverseLinearOperator(LinearOperator):
     """
-    Iterative solver for square linear system Ax=b, where x and b belong to (normed)
+    Abstract base class for the (approximate) inverse A_inv := A^{-1} of a
+    square matrix A. The result of A_inv.dot(b) is the (approximate) solution x
+    of the linear system A x = b, where x and b belong to the same (normed)
     vector space V.
 
+    We assume that the linear system is solved by an iterative method, which
+    needs a first guess `x0` and an exit condition based on `tol` and `maxiter`.
+
+    Concrete subclasses of this class must implement the `dot` method and take
+    care of any internal storage which might be necessary.
+
+    Parameters
+    ----------
+    A : psydac.linalg.basic.LinearOperator
+        Left-hand-side matrix A of linear system.
+        
+    x0 : psydac.linalg.basic.Vector
+        First guess of solution for iterative solver (optional).
+        
+    tol : float
+        Absolute tolerance for L2-norm of residual r = A*x - b.
+        
+    maxiter: int
+        Maximum number of iterations.
+        
+    verbose : bool
+        If True, L2-norm of residual r is printed at each iteration.
     """
 
-    @property
-    def space(self):
-        return self._space
+    def __init__(self, A, **kwargs):
+
+        assert isinstance(A, LinearOperator)
+        assert A.domain.dimension == A.codomain.dimension
+        domain = A.codomain
+        codomain = A.domain
+
+        if kwargs['x0'] is None:
+            kwargs['x0'] = codomain.zeros()
+
+        self._A = A
+        self._domain = domain
+        self._codomain = codomain
+
+        self._check_options(**kwargs)
+        self._options = kwargs
 
     @property
     def domain(self):
@@ -899,11 +949,40 @@ def dtype(self):
 
     @property
     def linop(self):
+        """
+        The linear operator L of which this object is the inverse L^{-1}.
+
+        The linear operator L can be modified in place, or replaced entirely
+        through the setter. A substitution should only be made in cases where
+        no other options are viable, as it breaks the one-to-one map between
+        the original linear operator L (passed to the constructor) and the
+        current `InverseLinearOperator` object L^{-1}. Use with extreme care!
+
+        """
         return self._A
+    
+    @linop.setter
+    def linop(self, a):
+        assert isinstance(a, LinearOperator)
+        assert a.domain is self.domain
+        assert a.codomain is self.codomain
+        self._A = a
 
-    @property
-    def options(self):
-        return self._options
+    def _check_options(self, **kwargs):
+        for key, value in kwargs.items():
+
+            if key == 'x0':
+                if value is not None:
+                    assert isinstance(value, Vector), "x0 must be a Vector or None"
+                    assert value.space == self.codomain, "x0 belongs to the wrong VectorSpace"
+            elif key == 'tol':
+                assert is_real(value), "tol must be a real number"
+                assert value > 0, "tol must be positive"
+            elif key == 'maxiter':
+                assert isinstance(value, int), "maxiter must be an int"
+                assert value > 0, "maxiter must be positive"
+            elif key == 'verbose':
+                assert isinstance(value, bool), "verbose must be a bool"
 
     def toarray(self):
         raise NotImplementedError('toarray() is not defined for InverseLinearOperators.')
@@ -914,76 +993,38 @@ def tosparse(self):
     def get_info(self):
         return self._info
 
-    def get_options(self):
-        return self._options.copy()
-
-    def set_options(self, **kwargs):
-        self._check_options(**kwargs)
-        self._options.update(kwargs)
+    def get_options(self, key=None):
+        """Get a copy of all the solver options, or a specific value of interest.
 
-    @abstractmethod
-    def _check_options(self, **kwargs):
-        pass
-
-    @abstractmethod
-    def transpose(self, conjugate=False):
-        pass
-
-    @staticmethod
-    def jacobi(A, b, out=None):
-        """
-        Jacobi preconditioner.
-
-        A : psydac.linalg.stencil.StencilMatrix | psydac.linalg.block.BlockLinearOperator
-            Left-hand-side matrix A of linear system.
-
-        b : psydac.linalg.stencil.StencilVector | psydac.linalg.block.BlockVector
-            Right-hand-side vector of linear system.
+        Parameters
+        ----------
+        key : str | None
+            Name of the specific option of interest (default: None).
 
         Returns
         -------
-        x : psydac.linalg.stencil.StencilVector | psydac.linalg.block.BlockVector
-            Preconditioner solution
+        dict | type(self._options['key']) | None
+            If `key` is given, get the specific option of interest. If there is
+            no such option, `None` is returned instead. If `key` is not given,
+            get a copy of all the solver options in a dictionary.
 
         """
-        from psydac.linalg.block   import BlockLinearOperator, BlockVector
-        from psydac.linalg.stencil import StencilMatrix, StencilVector
-
-        # In case A is None we return a zero vector
-        if A is None:
-            return b.space.zeros()
-
-        # Sanity checks
-        assert isinstance(A, (StencilMatrix, BlockLinearOperator))
-        assert isinstance(b, (StencilVector, BlockVector))
-        assert A.codomain.dimension == A.domain.dimension
-        assert A.codomain == b.space
-
-        #-------------------------------------------------------------
-        # Handle the case of a block linear system
-        if isinstance(A, BlockLinearOperator):
-            if out is not None:
-                for i, bi in enumerate(b.blocks):
-                    InverseLinearOperator.jacobi(A[i,i], bi, out=out[i])
-                return out
-            else:
-                x = [InverseLinearOperator.jacobi(A[i, i], bi) for i, bi in enumerate(b.blocks)]
-                y = BlockVector(b.space, blocks=x)
-                return y
-        #-------------------------------------------------------------
+        if key is None:
+            return self._options.copy()
+        else:
+            return self._options.get(key)
 
-        V = b.space
-        i = tuple(slice(s, e + 1) for s, e in zip(V.starts, V.ends))
+    def set_options(self, **kwargs):
+        """Set the solver options by passing keyword arguments.
+        """
+        self._check_options(**kwargs)
+        self._options.update(kwargs)
 
-        if out is not None:
-            b.copy(out=out)
-            out[i] /= A.diagonal()
-            out.update_ghost_regions()
-        else:
-            out = b.copy()
-            out[i] /= A.diagonal()
-            out.update_ghost_regions()
-            return out
+    def transpose(self, conjugate=False):
+        cls     = type(self)
+        At      = self.linop.transpose(conjugate=conjugate)
+        options = self._options
+        return cls(At, **options)
 
 #===============================================================================
 class LinearSolver(ABC):
diff --git a/psydac/linalg/block.py b/psydac/linalg/block.py
index eed6c7906..e8a9c44d6 100644
--- a/psydac/linalg/block.py
+++ b/psydac/linalg/block.py
@@ -777,6 +777,49 @@ def __sub__(self, M):
     #--------------------------------------
     # New properties/methods
     #--------------------------------------
+    def diagonal(self, *, inverse = False, out = None):
+        """Get the coefficients on the main diagonal as another BlockLinearOperator object.
+
+        Parameters
+        ----------
+        inverse : bool
+            If True, get the inverse of the diagonal. (Default: False).
+
+        out : BlockLinearOperator
+            If provided, write the diagonal entries into this matrix. (Default: None).
+
+        Returns
+        -------
+        BlockLinearOperator
+            The matrix which contains the main diagonal of self (or its inverse).
+
+        """
+        # Determine domain and codomain of result
+        V, W = self.domain, self.codomain
+        if inverse:
+            V, W = W, V
+
+        # Check the `out` argument, if `None` create a new BlockLinearOperator
+        if out is not None:
+            assert isinstance(out, BlockLinearOperator)
+            assert out.domain is V
+            assert out.codomain is W
+
+            # Set any off-diagonal blocks to zero
+            for i, j in out.nonzero_block_indices:
+                if i != j:
+                    out[i, j] = None
+        else:
+            out = BlockLinearOperator(V, W)
+
+        # Store the diagonal (or its inverse) into `out`
+        for i, j in self.nonzero_block_indices:
+            if i == j:
+                out[i, i] = self[i, i].diagonal(inverse = inverse, out = out[i, i])
+
+        return out
+
+    # ...
     @property
     def blocks(self):
         """ Immutable 2D view (tuple of tuples) of the linear operator,
@@ -868,7 +911,8 @@ def __setitem__(self, key, value):
             assert value.codomain is self.codomain[i]
 
         self._blocks[i,j] = value
-    
+
+    # ...
     def transform(self, operation):
         """
         Applies an operation on each block in this BlockLinearOperator.
@@ -881,6 +925,7 @@ def transform(self, operation):
         blocks = {ij: operation(Bij) for ij, Bij in self._blocks.items()}
         return BlockLinearOperator(self.domain, self.codomain, blocks=blocks)
 
+    # ...
     def backend(self):
         return self._backend
 
@@ -909,7 +954,6 @@ def copy(self, out=None):
 
         out.set_backend(self._backend)
         
-
     # ...
     def __imul__(self, a):
         for Bij in self._blocks.values():
diff --git a/psydac/linalg/solvers.py b/psydac/linalg/solvers.py
index 93817ac51..aa0b5e21b 100644
--- a/psydac/linalg/solvers.py
+++ b/psydac/linalg/solvers.py
@@ -1,20 +1,28 @@
 # coding: utf-8
 """
-This module provides iterative solvers and precondionners.
+This module provides iterative solvers and preconditioners.
 
 """
-from math import sqrt
 import numpy as np
+from math import sqrt
 
-from psydac.linalg.basic     import Vector, LinearOperator, InverseLinearOperator, IdentityOperator, ScaledLinearOperator
+from psydac.utilities.utils  import is_real
 from psydac.linalg.utilities import _sym_ortho
+from psydac.linalg.basic     import (Vector, LinearOperator,
+        InverseLinearOperator, IdentityOperator, ScaledLinearOperator)
+
+__all__ = (
+    'ConjugateGradient',
+    'PConjugateGradient',
+    'BiConjugateGradient',
+    'BiConjugateGradientStabilized',
+    'PBiConjugateGradientStabilized',
+    'MinimumResidual',
+    'LSMR',
+    'GMRES'
+)
 
-__all__ = ('ConjugateGradient', 'PConjugateGradient', 'BiConjugateGradient', 'BiConjugateGradientStabilized', 'MinimumResidual', 'LSMR', 'GMRES')
-
-def is_real(x):
-    from numbers import Number
-    return isinstance(x, Number) and np.isrealobj(x) and not isinstance(x, bool)
-
+#===============================================================================
 def inverse(A, solver, **kwargs):
     """
     A function to create objects of all InverseLinearOperator subclasses.
@@ -34,21 +42,35 @@ def inverse(A, solver, **kwargs):
 
     solver : str
         Preferred iterative solver. Options are: 'cg', 'pcg', 'bicg',
-        'bicgstab', 'minres', 'lsmr', 'gmres'.
+        'bicgstab', 'pbicgstab', 'minres', 'lsmr', 'gmres'.
 
     Returns
     -------
     obj : psydac.linalg.basic.InverseLinearOperator
-        More specifically: Returns the chosen subclass, for example psydac.linalg.solvers.ConjugateGradient
-        A linear operator acting as the inverse of A.
+        A linear operator acting as the inverse of A, of the chosen subclass
+        (for example psydac.linalg.solvers.ConjugateGradient).
 
     """
+
+    # Map each possible value of the `solver` string with a specific
+    # `InverseLinearOperator` subclass in this module:
+    solvers_dict = {
+        'cg'       : ConjugateGradient,
+        'pcg'      : PConjugateGradient,
+        'bicg'     : BiConjugateGradient,
+        'bicgstab' : BiConjugateGradientStabilized,
+        'pbicgstab': PBiConjugateGradientStabilized,
+        'minres'   : MinimumResidual,
+        'lsmr'     : LSMR,
+        'gmres'    : GMRES,
+    }
+
     # Check solver input
-    solvers = ('cg', 'pcg', 'bicg', 'bicgstab', 'minres', 'lsmr', 'gmres')
-    if solver not in solvers:
+    if solver not in solvers_dict:
         raise ValueError(f"Required solver '{solver}' not understood.")
 
     assert isinstance(A, LinearOperator)
+
     if isinstance(A, IdentityOperator):
         return A
     elif isinstance(A, ScaledLinearOperator):
@@ -57,20 +79,9 @@ def inverse(A, solver, **kwargs):
         return A.linop
 
     # Instantiate object of correct solver class
-    if solver == 'cg':
-        obj = ConjugateGradient(A, **kwargs)
-    elif solver == 'pcg':
-        obj = PConjugateGradient(A, **kwargs)
-    elif solver == 'bicg':
-        obj = BiConjugateGradient(A, **kwargs)
-    elif solver == 'bicgstab':
-        obj = BiConjugateGradientStabilized(A, **kwargs)
-    elif solver == 'minres':
-        obj = MinimumResidual(A, **kwargs)
-    elif solver == 'lsmr':
-        obj = LSMR(A, **kwargs)
-    elif solver == 'gmres':
-        obj = GMRES(A, **kwargs)
+    cls = solvers_dict[solver]
+    obj = cls(A, **kwargs)
+
     return obj
 
 #===============================================================================
@@ -111,52 +122,13 @@ class ConjugateGradient(InverseLinearOperator):
     """
     def __init__(self, A, *, x0=None, tol=1e-6, maxiter=1000, verbose=False, recycle=False):
 
-        assert isinstance(A, LinearOperator)
-        assert A.domain.dimension == A.codomain.dimension
-        domain = A.codomain
-        codomain = A.domain
-
-        if x0 is not None:
-            assert isinstance(x0, Vector)
-            assert x0.space is codomain
-        else:
-            x0 = codomain.zeros()
-
-        self._A = A
-        self._domain = domain
-        self._codomain = codomain
-        self._solver = 'cg'
         self._options = {"x0":x0, "tol":tol, "maxiter":maxiter, "verbose":verbose, "recycle":recycle}
-        self._check_options(**self._options)
-        self._tmps = {key: domain.zeros() for key in ("v", "r", "p")}
+        
+        super().__init__(A, **self._options)
+        
+        self._tmps = {key: self.domain.zeros() for key in ("v", "r", "p")}
         self._info = None
 
-    def _check_options(self, **kwargs):
-        for key, value in kwargs.items():
-
-            if key == 'x0':
-                if value is not None:
-                    assert isinstance(value, Vector), "x0 must be a Vector or None"
-                    assert value.space == self._codomain, "x0 belongs to the wrong VectorSpace"
-            elif key == 'tol':
-                assert is_real(value), "tol must be a real number"
-                assert value > 0, "tol must be positive"
-            elif key == 'maxiter':
-                assert isinstance(value, int), "maxiter must be an int"
-                assert value > 0, "maxiter must be positive"
-            elif key == 'verbose':
-                assert isinstance(value, bool), "verbose must be a bool"
-            elif key == 'recycle':
-                assert isinstance(value, bool), "recycle must be a bool"
-            else:
-                raise ValueError(f"Key '{key}' not understood. See self._options for allowed keys.")
-
-    def transpose(self, conjugate=False):
-        At = self._A.transpose(conjugate=conjugate)
-        solver = self._solver
-        options = self._options
-        return inverse(At, solver, **options)
-
     def solve(self, b, out=None):
         """
         Conjugate gradient algorithm for solving linear system Ax=b.
@@ -267,97 +239,63 @@ class PConjugateGradient(InverseLinearOperator):
     A LinearOperator subclass. Objects of this class are meant to be created using :func:~`solvers.inverse`.
     The .dot (and also the .solve) function are based on a preconditioned conjugate gradient method.
 
-    Preconditioned Conjugate Gradient (PCG) solves the symetric positive definte
-    system Ax = b. It assumes that pc(r) returns the solution to Ps = r,
-    where P is positive definite.
+    The Preconditioned Conjugate Gradient (PCG) algorithm solves the linear
+    system A x = b where A is a symmetric and positive-definite matrix, i.e.
+    A = A^T and y A y > 0 for any vector y. The preconditioner P is a matrix
+    which approximates the inverse of A. The algorithm assumes that P is also
+    symmetric and positive definite.
+
+    Since this is a matrix-free iterative method, both A and P are provided as
+    `LinearOperator` objects which must implement the `dot` method.
 
     Parameters
     ----------
-    A : psydac.linalg.stencil.StencilMatrix
-        Left-hand-side matrix A of linear system
-
-    pc: str
-        Preconditioner for A, it should approximate the inverse of A.
-        Can currently only be:
-        * The strings 'jacobi' or 'weighted_jacobi'. (rather obsolete, supply a callable instead, if possible)(14.02.: test weighted_jacobi)
-        The following should also be possible
-        * None, i.e. not pre-conditioning (this calls the standard `cg` method)
-        * A LinearSolver object (in which case the out parameter is used)
-        * A callable with two parameters (A, r), where A is the LinearOperator from above, and r is the residual.
+    A : psydac.linalg.basic.LinearOperator
+        Left-hand-side matrix A of the linear system. This should be symmetric
+        and positive definite.
+
+    pc: psydac.linalg.basic.LinearOperator
+        Preconditioner which should approximate the inverse of A (optional).
+        Like A, the preconditioner should be symmetric and positive definite.
 
     x0 : psydac.linalg.basic.Vector
         First guess of solution for iterative solver (optional).
 
     tol : float
-        Absolute tolerance for L2-norm of residual r = A*x - b.
+        Absolute tolerance for L2-norm of residual r = A x - b. (Default: 1e-6)
 
     maxiter: int
-        Maximum number of iterations.
+        Maximum number of iterations. (Default: 1000)
 
     verbose : bool
-        If True, L2-norm of residual r is printed at each iteration.
+        If True, the L2-norm of the residual r is printed at each iteration.
+        (Default: False)
 
     recycle : bool
-        Stores a copy of the output in x0 to speed up consecutive calculations of slightly altered linear systems
+        If True, a copy of the output is stored in x0 to speed up consecutive
+        calculations of slightly altered linear systems. (Default: False)
 
     """
-    def __init__(self, A, *, pc='jacobi', x0=None, tol=1e-6, maxiter=1000, verbose=False, recycle=False):
+    def __init__(self, A, *, pc=None, x0=None, tol=1e-6, maxiter=1000, verbose=False, recycle=False):
 
-        assert isinstance(A, LinearOperator)
-        assert A.domain.dimension == A.codomain.dimension
-        domain = A.codomain
-        codomain = A.domain
-
-        if x0 is not None:
-            assert isinstance(x0, Vector)
-            assert x0.space is codomain
-        else:
-            x0 = codomain.zeros()
-
-        self._A = A
-        self._domain = domain
-        self._codomain = codomain
-        self._solver = 'pcg'
         self._options = {"x0":x0, "pc":pc, "tol":tol, "maxiter":maxiter, "verbose":verbose, "recycle":recycle}
-        self._check_options(**self._options)
-        tmps_codomain = {key: codomain.zeros() for key in ("p", "s")}
-        tmps_domain = {key: domain.zeros() for key in ("v", "r")}
+        
+        super().__init__(A, **self._options)
+        
+        if pc is None:
+            self._options['pc'] = IdentityOperator(self.domain)
+        else:
+            assert isinstance(pc, LinearOperator)
+            
+        tmps_codomain = {key: self.codomain.zeros() for key in ("p", "s")}
+        tmps_domain = {key: self.domain.zeros() for key in ("v", "r")}
         self._tmps = {**tmps_codomain, **tmps_domain}
         self._info = None
 
-    def _check_options(self, **kwargs):
-        for key, value in kwargs.items():
-
-            if key == 'pc':
-                assert value is not None, "pc may not be None"
-                assert value == 'jacobi', "unsupported preconditioner"
-            elif key == 'x0':
-                if value is not None:
-                    assert isinstance(value, Vector), "x0 must be a Vector or None"
-                    assert value.space == self._codomain, "x0 belongs to the wrong VectorSpace"
-            elif key == 'tol':
-                assert is_real(value), "tol must be a real number"
-                assert value > 0, "tol must be positive"
-            elif key == 'maxiter':
-                assert isinstance(value, int), "maxiter must be an int"
-                assert value > 0, "maxiter must be positive"
-            elif key == 'verbose':
-                assert isinstance(value, bool), "verbose must be a bool"
-            elif key == 'recycle':
-                assert isinstance(value, bool), "recycle must be a bool"
-            else:
-                raise ValueError(f"Key '{key}' not understood. See self._options for allowed keys.")
-
-    def transpose(self, conjugate=False):
-        At = self._A.transpose(conjugate=conjugate)
-        solver = self._solver
-        options = self._options
-        return inverse(At, solver, **options)
-
     def solve(self, b, out=None):
         """
         Preconditioned Conjugate Gradient (PCG) solves the symetric positive definte
-        system Ax = b. It assumes that pc(r) returns the solution to Ps = r,
+        system Ax = b. It assumes that pc.dot(r) returns the solution to Ps = r,
         where P is positive definite.
         Info can be accessed using get_info(), see :func:~`basic.InverseLinearOperator.get_info`.
 
@@ -390,6 +328,8 @@ def solve(self, b, out=None):
 
         assert isinstance(b, Vector)
         assert b.space is domain
+    
+        assert isinstance(pc, LinearOperator)
 
         # First guess of solution
         if out is not None:
@@ -398,22 +338,6 @@ def solve(self, b, out=None):
 
         x = x0.copy(out=out)
 
-        # Preconditioner
-        assert pc is not None
-        if pc == 'jacobi':
-            psolve = lambda r, out: InverseLinearOperator.jacobi(A, r, out)
-        #elif pc == 'weighted_jacobi':
-        #    psolve = lambda r, out: InverseLinearOperator.weighted_jacobi(A, r, out) # allows for further specification not callable like this!
-        #elif isinstance(pc, str):
-        #    pcfun = getattr(InverseLinearOperator, str)
-        #    #pcfun = globals()[pc]
-        #    psolve = lambda r: pcfun(A, r)
-        #elif isinstance(pc, LinearSolver):
-        #    s = b.space.zeros()
-        #    psolve = lambda r: pc.solve(r, out=s)
-        #elif hasattr(pc, '__call__'):
-        #    psolve = lambda r: pc(A, r)
-
         # Extract local storage
         v = self._tmps["v"]
         r = self._tmps["r"]
@@ -425,7 +349,7 @@ def solve(self, b, out=None):
         b.copy(out=r)
         r       -= v
         nrmr_sqr = r.dot(r).real
-        psolve(r, out=s)
+        pc.dot(r, out=s)
         am       = s.dot(r)
         s.copy(out=p)
 
@@ -453,7 +377,7 @@ def solve(self, b, out=None):
             r.mul_iadd(-l, v) # this is r -= l*v
 
             nrmr_sqr = r.dot(r).real
-            psolve(r, out=s)
+            pc.dot(r, out=s)
 
             am1 = s.dot(r)
 
@@ -518,53 +442,14 @@ class BiConjugateGradient(InverseLinearOperator):
     """
     def __init__(self, A, *, x0=None, tol=1e-6, maxiter=1000, verbose=False, recycle=False):
 
-        assert isinstance(A, LinearOperator)
-        assert A.domain.dimension == A.codomain.dimension
-        domain = A.codomain
-        codomain = A.domain
-
-        if x0 is not None:
-            assert isinstance(x0, Vector)
-            assert x0.space is codomain
-        else:
-            x0 = codomain.zeros()
-
-        self._A = A
-        self._Ah = A.H
-        self._domain = domain
-        self._codomain = codomain
-        self._solver = 'bicg'
         self._options = {"x0":x0, "tol":tol, "maxiter":maxiter, "verbose":verbose, "recycle":recycle}
-        self._check_options(**self._options)
-        self._tmps = {key: domain.zeros() for key in ("v", "r", "p", "vs", "rs", "ps")}
+        
+        super().__init__(A, **self._options)
+        
+        self._Ah = A.H
+        self._tmps = {key: self.domain.zeros() for key in ("v", "r", "p", "vs", "rs", "ps")}
         self._info = None
 
-    def _check_options(self, **kwargs):
-        for key, value in kwargs.items():
-
-            if key == 'x0':
-                if value is not None:
-                    assert isinstance(value, Vector), "x0 must be a Vector or None"
-                    assert value.space == self._codomain, "x0 belongs to the wrong VectorSpace"
-            elif key == 'tol':
-                assert is_real(value), "tol must be a real number"
-                assert value > 0, "tol must be positive"
-            elif key == 'maxiter':
-                assert isinstance(value, int), "maxiter must be an int"
-                assert value > 0, "maxiter must be positive"
-            elif key == 'verbose':
-                assert isinstance(value, bool), "verbose must be a bool"
-            elif key == 'recycle':
-                assert isinstance(value, bool), "recycle must be a bool"
-            else:
-                raise ValueError(f"Key '{key}' not understood. See self._options for allowed keys.")
-
-    def transpose(self, conjugate=False):
-        At = self._A.transpose(conjugate=conjugate)
-        solver = self._solver
-        options = self._options
-        return inverse(At, solver, **options)
-
     def solve(self, b, out=None):
         """
         Biconjugate gradient (BCG) algorithm for solving linear system Ax=b.
@@ -746,52 +631,13 @@ class BiConjugateGradientStabilized(InverseLinearOperator):
     """
     def __init__(self, A, *, x0=None, tol=1e-6, maxiter=1000, verbose=False, recycle=False):
 
-        assert isinstance(A, LinearOperator)
-        assert A.domain.dimension == A.codomain.dimension
-        domain = A.codomain
-        codomain = A.domain
-
-        if x0 is not None:
-            assert isinstance(x0, Vector)
-            assert x0.space is codomain
-        else:
-            x0 = codomain.zeros()
-
-        self._A = A
-        self._domain = domain
-        self._codomain = codomain
-        self._solver = 'bicgstab'
         self._options = {"x0": x0, "tol": tol, "maxiter": maxiter, "verbose": verbose, "recycle":recycle}
-        self._check_options(**self._options)
-        self._tmps = {key: domain.zeros() for key in ("v", "r", "p", "vr", "r0")}
+        
+        super().__init__(A, **self._options)
+        
+        self._tmps = {key: self.domain.zeros() for key in ("v", "r", "p", "vr", "r0")}
         self._info = None
 
-    def _check_options(self, **kwargs):
-        for key, value in kwargs.items():
-
-            if key == 'x0':
-                if value is not None:
-                    assert isinstance(value, Vector), "x0 must be a Vector or None"
-                    assert value.space == self._codomain, "x0 belongs to the wrong VectorSpace"
-            elif key == 'tol':
-                assert is_real(value), "tol must be a real number"
-                assert value > 0, "tol must be positive"
-            elif key == 'maxiter':
-                assert isinstance(value, int), "maxiter must be an int"
-                assert value > 0, "maxiter must be positive"
-            elif key == 'verbose':
-                assert isinstance(value, bool), "verbose must be a bool"
-            elif key == 'recycle':
-                assert isinstance(value, bool), "recycle must be a bool"
-            else:
-                raise ValueError(f"Key '{key}' not understood. See self._options for allowed keys.")
-
-    def transpose(self, conjugate=False):
-        At = self._A.transpose(conjugate=conjugate)
-        solver = self._solver
-        options = self._options
-        return inverse(At, solver, **options)
-
     def solve(self, b, out=None):
         """
         Biconjugate gradient stabilized method (BCGSTAB) algorithm for solving linear system Ax=b.
@@ -944,6 +790,232 @@ def solve(self, b, out=None):
     def dot(self, b, out=None):
         return self.solve(b, out=out)
 
+#===============================================================================
+class PBiConjugateGradientStabilized(InverseLinearOperator):
+    """
+    A LinearOperator subclass. Objects of this class are meant to be created using :func:~`solvers.inverse`.
+    The .dot (and also the .solve) function are based on the
+    preconditioned Biconjugate gradient Stabilized (PBCGSTAB) algorithm for solving linear system Ax=b.
+    Implementation from [1], page 251.
+    Parameters
+    ----------
+    A : psydac.linalg.basic.LinearOperator
+        Left-hand-side matrix A of linear system; individual entries A[i,j]
+        can't be accessed, but A has 'shape' attribute and provides 'dot(p)'
+        function (i.e. matrix-vector product A*p).
+    pc: psydac.linalg.basic.LinearOperator
+        Preconditioner for A, it should approximate the inverse of A (can be None).
+    x0 : psydac.linalg.basic.Vector
+        First guess of solution for iterative solver (optional).
+    tol : float
+        Absolute tolerance for 2-norm of residual r = A*x - b.
+    maxiter: int
+        Maximum number of iterations.
+    verbose : bool
+        If True, 2-norm of residual r is printed at each iteration.
+    References
+    ----------
+    [1] A. Maister, Numerik linearer Gleichungssysteme, Springer ed. 2015.
+    """
+    def __init__(self, A, *, pc=None, x0=None, tol=1e-6, maxiter=1000, verbose=False, recycle=False):
+
+        self._options = {"pc": pc, "x0": x0, "tol": tol, "maxiter": maxiter, "verbose": verbose, "recycle": recycle}
+
+        super().__init__(A, **self._options)
+
+        if pc is None:
+            self._options['pc'] = IdentityOperator(self.domain)
+        else:
+            assert isinstance(pc, LinearOperator)
+
+        self._tmps = {key: self.domain.zeros() for key in ("v", "r", "s", "t", 
+                                                      "vp", "rp", "sp", "tp",
+                                                      "pp", "av", "app", "osp", 
+                                                      "rp0")}
+        self._info = None
+
+    def solve(self, b, out=None):
+        """
+        Preconditioned biconjugate gradient stabilized method (PBCGSTAB) algorithm for solving linear system Ax=b.
+        Implementation from [1], page 251.
+        Parameters
+        ----------
+        b : psydac.linalg.basic.Vector
+            Right-hand-side vector of linear system. Individual entries b[i] need
+            not be accessed, but b has 'shape' attribute and provides 'copy()' and
+            'dot(p)' functions (dot(p) is the vector inner product b*p ); moreover,
+            scalar multiplication and sum operations are available.
+        out : psydac.linalg.basic.Vector | NoneType
+            The output vector, or None (optional).
+        Returns
+        -------
+        x : psydac.linalg.basic.Vector
+            Numerical solution of linear system. To check the convergence of the solver,
+            use the method InverseLinearOperator.get_info().
+        
+        info : dict
+            Dictionary containing convergence information:
+              - 'niter'    = (int) number of iterations
+              - 'success'  = (boolean) whether convergence criteria have been met
+              - 'res_norm' = (float) 2-norm of residual vector r = A*x - b.
+        References
+        ----------
+        [1] A. Maister, Numerik linearer Gleichungssysteme, Springer ed. 2015.
+        """
+
+        A = self._A
+        domain = self._domain
+        codomain = self._codomain
+        options = self._options
+        pc = options["pc"]
+        x0 = options["x0"]
+        tol = options["tol"]
+        maxiter = options["maxiter"]
+        verbose = options["verbose"]
+        recycle = options["recycle"]
+
+        assert isinstance(b, Vector)
+        assert b.space is domain
+
+        assert isinstance(pc, LinearOperator)
+
+        # first guess of solution
+        if out is not None:
+            assert isinstance(out, Vector)
+            assert out.space == codomain
+            out *= 0
+            if x0 is None:
+                x = out
+            else:
+                assert x0.shape == (A.shape[0],)
+                out += x0
+                x = out
+        else:
+            if x0 is None:
+                x = b.copy()
+                x *= 0.0
+            else:
+                assert x0.shape == (A.shape[0],)
+                x = x0.copy()
+
+        # preconditioner (must have a .solve method)
+        assert isinstance(pc, LinearOperator)
+
+        # extract temporary vectors
+        v = self._tmps['v']
+        r = self._tmps['r']
+        s = self._tmps['s']
+        t = self._tmps['t']
+
+        vp = self._tmps['vp']
+        rp = self._tmps['rp']
+        pp = self._tmps['pp']
+        sp = self._tmps['sp']
+        tp = self._tmps['tp']
+
+        av = self._tmps['av']
+
+        app = self._tmps['app']
+        osp = self._tmps['osp']
+
+        # first values: r = b - A @ x, rp = pp = PC @ r, rhop = |rp|^2
+        A.dot(x, out=v)
+        b.copy(out=r)
+        r -= v
+
+        pc.dot(r, out=rp)
+        rp.copy(out=pp)
+
+        rhop = rp.dot(rp)
+
+        # save initial residual vector rp0
+        rp0 = self._tmps['rp0']
+        rp.copy(out=rp0)
+
+        # squared residual norm and squared tolerance
+        res_sqr = r.dot(r)
+        tol_sqr = tol**2
+
+        if verbose:
+            print("Pre-conditioned BICGSTAB solver:")
+            print("+---------+---------------------+")
+            print("+ Iter. # | L2-norm of residual |")
+            print("+---------+---------------------+")
+            template = "| {:7d} | {:19.2e} |"
+
+        # iterate to convergence or maximum number of iterations
+        niter = 0
+
+        while res_sqr > tol_sqr and niter < maxiter:
+
+            # v = A @ pp, vp = PC @ v, alphap = rhop/(vp.rp0)
+            A.dot(pp, out=v)
+            pc.dot(v, out=vp)
+            alphap = rhop / vp.dot(rp0)
+
+            # s = r - alphap*v, sp = PC @ s
+            r.copy(out=s)
+            v.copy(out=av)
+            av *= alphap
+            s -= av
+            pc.dot(s, out=sp)
+
+            # t = A @ sp, tp = PC @ t, omegap = (tp.sp)/(tp.tp)
+            A.dot(sp, out=t)
+            pc.dot(t, out=tp)
+            omegap = tp.dot(sp) / tp.dot(tp)
+
+            # x = x + alphap*pp + omegap*sp
+            pp.copy(out=app)
+            sp.copy(out=osp)
+            app *= alphap
+            osp *= omegap
+            x += app
+            x += osp
+
+            # r = s - omegap*t, rp = sp - omegap*tp
+            s.copy(out=r)
+            t *= omegap
+            r -= t
+
+            sp.copy(out=rp)
+            tp *= omegap
+            rp -= tp
+
+            # rhop_new = rp.rp0, betap = (alphap*rhop_new)/(omegap*rhop)
+            rhop_new = rp.dot(rp0)
+            betap = (alphap*rhop_new) / (omegap*rhop)
+            rhop = 1*rhop_new
+
+            # pp = rp + betap*(pp - omegap*vp)
+            vp *= omegap
+            pp -= vp
+            pp *= betap
+            pp += rp
+
+            # new residual norm
+            res_sqr = r.dot(r)
+
+            niter += 1
+
+            if verbose:
+                print(template.format(niter, sqrt(res_sqr)))
+
+        if verbose:
+            print("+---------+---------------------+")
+
+        # convergence information
+        self._info = {'niter': niter, 'success': res_sqr <
+                tol_sqr, 'res_norm': sqrt(res_sqr)}
+
+        if recycle:
+            x.copy(out=self._options["x0"])
+
+        return x
+
+    def dot(self, b, out=None):
+        return self.solve(b, out=out)
+
 #===============================================================================
 class MinimumResidual(InverseLinearOperator):
     """
@@ -992,52 +1064,12 @@ class MinimumResidual(InverseLinearOperator):
     """
     def __init__(self, A, *, x0=None, tol=1e-6, maxiter=1000, verbose=False, recycle=False):
 
-        assert isinstance(A, LinearOperator)
-        assert A.domain.dimension == A.codomain.dimension
-        assert A.domain.dtype == float
-        domain = A.codomain
-        codomain = A.domain
-
-        if x0 is not None:
-            assert isinstance(x0, Vector)
-            assert x0.space is codomain
-        else:
-            x0 = codomain.zeros()
-
-        self._A = A
-        self._domain = domain
-        self._codomain = codomain
-        self._solver = 'minres'
         self._options = {"x0":x0, "tol":tol, "maxiter":maxiter, "verbose":verbose, "recycle":recycle}
-        self._check_options(**self._options)
-        self._tmps = {key: domain.zeros() for key in ("res_old", "res_new", "w_new", "w_work", "w_old", "v", "y")}
-        self._info = None
-
-    def _check_options(self, **kwargs):
-        for key, value in kwargs.items():
-
-            if key == 'x0':
-                if value is not None:
-                    assert isinstance(value, Vector), "x0 must be a Vector or None"
-                    assert value.space == self._codomain, "x0 belongs to the wrong VectorSpace"
-            elif key == 'tol':
-                assert is_real(value), "tol must be a real number"
-                assert value > 0, "tol must be positive"
-            elif key == 'maxiter':
-                assert isinstance(value, int), "maxiter must be an int"
-                assert value > 0, "maxiter must be positive"
-            elif key == 'verbose':
-                assert isinstance(value, bool), "verbose must be a bool"
-            elif key == 'recycle':
-                assert isinstance(value, bool), "recycle must be a bool"
-            else:
-                raise ValueError(f"Key '{key}' not understood. See self._options for allowed keys.")
 
-    def transpose(self, conjugate=False):
-        At = self._A.transpose(conjugate=conjugate)
-        solver = self._solver
-        options = self._options
-        return inverse(At, solver, **options)
+        super().__init__(A, **self._options)
+        
+        self._tmps = {key: self.domain.zeros() for key in ("res_old", "res_new", "w_new", "w_work", "w_old", "v", "y")}
+        self._info = None
 
     def solve(self, b, out=None):
         """
@@ -1324,67 +1356,30 @@ class LSMR(InverseLinearOperator):
     """
     def __init__(self, A, *, x0=None, tol=None, atol=None, btol=None, maxiter=1000, conlim=1e8, verbose=False, recycle=False):
 
-        assert isinstance(A, LinearOperator)
-        assert A.domain.dimension == A.codomain.dimension
-        domain = A.codomain
-        codomain = A.domain
-
-        if x0 is not None:
-            assert isinstance(x0, Vector)
-            assert x0.space is codomain
-        else:
-            x0 = codomain.zeros()
-
-        self._A = A
-        self._domain = domain
-        self._codomain = codomain
-        self._solver = 'lsmr'
         self._options = {"x0":x0, "tol":tol, "atol":atol, "btol":btol,
                          "maxiter":maxiter, "conlim":conlim, "verbose":verbose, "recycle":recycle}
-        self._check_options(**self._options)
+        
+        super().__init__(A, **self._options)
+
+        # check additional options
+        if atol is not None:
+            assert is_real(atol), "atol must be a real number"
+            assert atol >= 0, "atol must not be negative"
+        if btol is not None:
+            assert is_real(btol), "btol must be a real number"
+            assert btol >= 0, "btol must not be negative"
+        assert is_real(conlim), "conlim must be a real number" # actually an integer?
+        assert conlim > 0, "conlim must be positive" # supposedly
+        
         self._info = None
         self._successful = None
-        tmps_domain = {key: domain.zeros() for key in ("u", "u_work")}
-        tmps_codomain = {key: codomain.zeros() for key in ("v", "v_work", "h", "hbar")}
+        tmps_domain = {key: self.domain.zeros() for key in ("u", "u_work")}
+        tmps_codomain = {key: self.codomain.zeros() for key in ("v", "v_work", "h", "hbar")}
         self._tmps = {**tmps_codomain, **tmps_domain}
 
     def get_success(self):
         return self._successful
 
-    def _check_options(self, **kwargs):
-        for key, value in kwargs.items():
-
-            if key == 'x0':
-                if value is not None:
-                    assert isinstance(value, Vector), "x0 must be a Vector or None"
-                    assert value.space == self._codomain, "x0 belongs to the wrong VectorSpace"
-            elif key == 'tol':
-                if value is not None:
-                    assert is_real(value), "tol must be a real number"
-                    assert value > 0, "tol must be positive" # suppose atol/btol must also be positive numbers
-            elif key == 'atol' or key == 'btol':
-                if value is not None:
-                    assert is_real(value), "atol/btol must be a real number"
-                    assert value >= 0, "atol/btol must not be negative"
-            elif key == 'maxiter':
-                assert isinstance(value, int), "maxiter must be an int"
-                assert value > 0, "maxiter must be positive"
-            elif key == 'conlim':
-                assert is_real(value), "conlim must be a real number" # actually an integer?
-                assert value > 0, "conlim must be positive" # supposedly
-            elif key == 'verbose':
-                assert isinstance(value, bool), "verbose must be a bool"
-            elif key == 'recycle':
-                assert isinstance(value, bool), "recycle must be a bool"
-            else:
-                raise ValueError(f"Key '{key}' not understood. See self._options for allowed keys.")
-
-    def transpose(self, conjugate=False):
-        At = self._A.transpose(conjugate=conjugate)
-        solver = self._solver
-        options = self._options
-        return inverse(At, solver, **options)
-
     def solve(self, b, out=None):
         """Iterative solver for least-squares problems.
         lsmr solves the system of linear equations ``Ax = b``. If the system
@@ -1711,56 +1706,17 @@ class GMRES(InverseLinearOperator):
     """
     def __init__(self, A, *, x0=None, tol=1e-6, maxiter=100, verbose=False, recycle=False):
 
-        assert isinstance(A, LinearOperator)
-        assert A.domain.dimension == A.codomain.dimension
-        domain = A.codomain
-        codomain = A.domain
+        self._options = {"x0":x0, "tol":tol, "maxiter":maxiter, "verbose":verbose, "recycle":recycle}
 
-        if x0 is not None:
-            assert isinstance(x0, Vector)
-            assert x0.space is codomain
-        else:
-            x0 = codomain.zeros()
+        super().__init__(A, **self._options)
 
-        self._A = A
-        self._domain = domain
-        self._codomain = codomain
-        self._solver = 'gmres'
-        self._options = {"x0":x0, "tol":tol, "maxiter":maxiter, "verbose":verbose, "recycle":recycle}
-        self._check_options(**self._options) 
-        self._tmps = {key: domain.zeros() for key in ("r", "p")}
+        self._tmps = {key: self.domain.zeros() for key in ("r", "p")}
 
         # Initialize upper Hessenberg matrix
         self._H = np.zeros((self._options["maxiter"] + 1, self._options["maxiter"]), dtype=A.domain.dtype)
         self._Q = []
         self._info = None
 
-    def _check_options(self, **kwargs):
-        for key, value in kwargs.items():
-
-            if key == 'x0':
-                if value is not None:
-                    assert isinstance(value, Vector), "x0 must be a Vector or None"
-                    assert value.space == self._codomain, "x0 belongs to the wrong VectorSpace"
-            elif key == 'tol':
-                assert is_real(value), "tol must be a real number"
-                assert value > 0, "tol must be positive"
-            elif key == 'maxiter':
-                assert isinstance(value, int), "maxiter must be an int"
-                assert value > 0, "maxiter must be positive"
-            elif key == 'verbose':
-                assert isinstance(value, bool), "verbose must be a bool"
-            elif key == 'recycle':
-                assert isinstance(value, bool), "recycle must be a bool"
-            else:
-                raise ValueError(f"Key '{key}' not understood. See self._options for allowed keys.")
-
-    def transpose(self, conjugate=False):
-        At = self._A.transpose(conjugate=conjugate)
-        solver = self._solver
-        options = self._options
-        return inverse(At, solver, **options)
-
     def solve(self, b, out=None):
         """
         Generalized minimal residual algorithm for solving linear system Ax=b.
@@ -1928,4 +1884,4 @@ def apply_givens_rotation(self, k, sn, cn):
 
     def dot(self, b, out=None):
         return self.solve(b, out=out)
-    
\ No newline at end of file
+
diff --git a/psydac/linalg/stencil.py b/psydac/linalg/stencil.py
index cf4fcddc1..2c8c6fa96 100644
--- a/psydac/linalg/stencil.py
+++ b/psydac/linalg/stencil.py
@@ -8,7 +8,7 @@
 import numpy as np
 
 from types        import MappingProxyType
-from scipy.sparse import coo_matrix
+from scipy.sparse import coo_matrix, diags as sp_diags
 from mpi4py       import MPI
 
 from psydac.linalg.basic  import VectorSpace, Vector, LinearOperator
@@ -1368,30 +1368,61 @@ def _exchange_assembly_data_serial(self):
 
                 # Copy data from left to right
                 idx_to   = tuple( idx_front + [slice( m*p, m*p+p)] + idx_back )
-                idx_from = tuple( idx_front + [ slice(-m*p,-m*p+p) if (-m*p+p)!=0 else slice(-m*p,None)] + idx_back )
+                idx_from = tuple( idx_front + [slice(-m*p,-m*p+p) if (-m*p+p)!=0 else slice(-m*p,None)] + idx_back )
                 self._data[idx_to] += self._data[idx_from]
 
-    def diagonal(self):
-        if self._diag_indices is None:
-            cm    = self.codomain.shifts
-            dm    = self.domain.shifts
-            ss    = self.codomain.starts
-            pp    = [compute_diag_len(p,mj,mi)-(p+1) for p,mi,mj in zip(self._pads, cm, dm)]
-            nrows = tuple(e-s+1 for s,e in zip(self.codomain.starts, self.codomain.ends))
-            indices = [np.zeros(np.product(nrows), dtype=int) for _ in range(2*len(nrows))]
-            l = 0
-            for xx in np.ndindex(*nrows):
-                ii = [m*p+x for m,p,x in zip(self.domain.shifts, self.domain.pads, xx)]
-                jj = [p+x+s-((x+s)//mi)*mj for (x,mi,mj,p,s) in zip(xx,cm,dm,pp,ss)]
-                for k in range(len(nrows)):
-                    indices[k][l] = ii[k]
-                    indices[k+len(nrows)][l] = jj[k]
-                l += 1
-            self._diag_indices = tuple(indices)
+    # ...
+    def diagonal(self, *, inverse = False, out = None):
+        """
+        Get the coefficients on the main diagonal as a StencilDiagonalMatrix object.
+
+        Parameters
+        ----------
+        inverse : bool
+            If True, get the inverse of the diagonal. (Default: False).
+
+        out : StencilDiagonalMatrix
+            If provided, write the diagonal entries into this matrix. (Default: None).
+
+        Returns
+        -------
+        StencilDiagonalMatrix
+            The matrix which contains the main diagonal of self (or its inverse).
+
+        """
+        # Check `inverse` argument
+        assert isinstance(inverse, bool)
+
+        # Determine domain and codomain of the StencilDiagonalMatrix
+        V, W = self.domain, self.codomain
+        if inverse:
+            V, W = W, V
+
+        # Check `out` argument
+        if out is not None:
+            assert isinstance(out, StencilDiagonalMatrix)
+            assert out.domain is V
+            assert out.codomain is W
+
+
+        # Extract diagonal data from self and identify output array
+        diagonal_indices = self._get_diagonal_indices()
+        diag = self._data[diagonal_indices]
+        data = out._data if out else None
+
+        # Calculate entries of StencilDiagonalMatrix
+        if inverse:
+            data = np.divide(1, diag, out=data)
+        elif out:
+            np.copyto(data, diag)
         else:
-            nrows   = tuple(e-s+1 for s,e in zip(self.codomain.starts, self.codomain.ends))
+            data = diag.copy()
+
+        # If needed create a new StencilDiagonalMatrix object
+        if out is None:
+            out = StencilDiagonalMatrix(V, W, data)
 
-        return self._data[self._diag_indices].reshape(nrows)
+        return out
 
     # ...
     def topetsc(self):
@@ -1813,6 +1844,224 @@ def set_backend(self, backend):
             self._args.pop('ndiags')
             self._func = dot.func
 
+    # ...
+    def _get_diagonal_indices(self):
+        """
+        Compute the indices which should be applied to self._data in order to
+        get the matrix entries on the main diagonal. The result is also stored
+        in self._diag_indices, and retrieved from there on successive calls.
+
+        Returns
+        -------
+        tuple[numpy.ndarray, ndim]
+            The diagonal indices as a tuple of NumPy arrays of identical shape
+            (n1, n2, n3, ...).
+
+        """
+
+        if self._diag_indices is None:
+
+            dp    = self.domain.pads
+            dm    = self.domain.shifts
+            cm    = self.codomain.shifts
+            ss    = self.codomain.starts
+            pp    = [compute_diag_len(p, mj, mi) - p - 1 for p, mi, mj in zip(self._pads, cm, dm)]
+            nrows = [e - s + 1 for s, e in zip(self.codomain.starts, self.codomain.ends)]
+            ndim  = self.domain.ndim
+
+            indices = [np.zeros(np.product(nrows), dtype=int) for _ in range(2 * ndim)]
+
+            for l, xx in enumerate(np.ndindex(*nrows)):
+                ii = [m * p + x for m, p, x in zip(dm, dp, xx)]
+                jj = [p + x + s - ((x+s) // mi) * mj for x, mi, mj, p, s in zip(xx, cm, dm, pp, ss)]
+                for k in range(ndim):
+                    indices[k][l] = ii[k]
+                    indices[k + ndim][l] = jj[k]
+
+            self._diag_indices = tuple(idx.reshape(nrows) for idx in indices)
+
+        return self._diag_indices
+
+#===============================================================================
+class StencilDiagonalMatrix(LinearOperator):
+    """
+    Linear operator which operates between stencil vector spaces, and which can
+    be represented by a matrix with non-zero entries only on its main diagonal.
+    As such this operator is completely local and requires no data communication.
+
+    We assume that the vectors in the domain and the codomain have the same
+    shape and are distributed in the same way.
+
+    Parameters
+    ----------
+    V : psydac.linalg.stencil.StencilVectorSpace
+        Domain of the new linear operator.
+
+    W : psydac.linalg.stencil.StencilVectorSpace
+        Codomain of the new linear operator.
+
+    """
+    def __init__(self, V, W, data):
+
+        # Check domain and codomain
+        assert isinstance(V, StencilVectorSpace)
+        assert isinstance(W, StencilVectorSpace)
+        assert V.starts == W.starts
+        assert V.ends   == W.ends
+
+        data = np.asarray(data)
+
+        # Check shape of provided data
+        shape = tuple(e - s + 1 for s, e in zip(V.starts, V.ends))
+        assert data.shape == shape
+
+        # Store info in object
+        self._domain   = V
+        self._codomain = W
+        self._data     = data
+
+    #--------------------------------------
+    # Abstract interface
+    #--------------------------------------
+    @property
+    def domain(self):
+        return self._domain
+
+    @property
+    def codomain(self):
+        return self._codomain
+
+    @property
+    def dtype(self):
+        return self._data.dtype
+
+    def tosparse(self):
+        return sp_diags(self._data.ravel())
+
+    def toarray(self):
+        return self._data.copy()
+
+    def dot(self, v, out=None):
+
+        assert isinstance(v, StencilVector)
+        assert v.space is self.domain
+
+        if out is not None:
+            assert isinstance(out, StencilVector)
+            assert out.space is self.codomain
+        else:
+            out = self.codomain.zeros()
+
+        V = self.domain
+        i = tuple(slice(s, e + 1) for s, e in zip(V.starts, V.ends))
+        np.multiply(self._data, v[i], out=out[i])
+
+        out.ghost_regions_in_sync = False
+
+        return out
+
+    # ...
+    # TODO [YG 22.01.2024]: idot function will require a dedicated kernel
+    # ...
+
+    def transpose(self, *, conjugate=False, out=None):
+
+        assert isinstance(conjugate, bool)
+
+        if out is not None:
+            assert isinstance(out, StencilDiagonalMatrix)
+            assert out.domain is self.codomain
+            assert out.codomain is self.domain
+
+        if not (conjugate and self.dtype is complex):
+
+            if out is None:
+                data = self._data.copy()
+            else:
+                np.copyto(out._data, self._data, casting='no')
+
+        else:
+
+            if out is None:
+                data = np.conjugate(self._data, casting='no')
+            else:
+                np.conjugate(self._data, out=out._data, casting='no')
+
+        if out is None:
+            out = StencilDiagonalMatrix(self.codomain, self.domain, data)
+
+        return out
+
+    #--------------------------------------
+    # Other properties/methods
+    #--------------------------------------
+    def copy(self, *, out=None):
+
+        if out is self:
+            return self
+
+        if out is None:
+            data = self._data.copy()
+            out = StencilDiagonalMatrix(self.domain, self.codomain, data)
+        else:
+            assert isinstance(out, StencilDiagonalMatrix)
+            assert out.domain is self.domain
+            assert out.codomain is self.codomain
+            np.copyto(out._data, self._data, casting='no')
+
+        return out
+
+    def diagonal(self, *, inverse = False, out = None):
+        """
+        Get the coefficients on the main diagonal as a StencilDiagonalMatrix object.
+
+        In the default case (inverse=False, out=None) self is returned.
+
+        Parameters
+        ----------
+        inverse : bool
+            If True, get the inverse of the diagonal. (Default: False).
+
+        out : StencilDiagonalMatrix
+            If provided, write the diagonal entries into this matrix. (Default: None).
+
+        Returns
+        -------
+        StencilDiagonalMatrix
+            Either self, or another StencilDiagonalMatrix with the diagonal inverse.
+
+        """
+        # Check `inverse` argument
+        assert isinstance(inverse, bool)
+
+        # Determine domain and codomain of the `out` matrix
+        V, W = self.domain, self.codomain
+        if inverse:
+            V, W = W, V
+
+        # Check `out` argument and identify `data` array of output vector
+        if out is None:
+            data = None
+        else:
+            assert isinstance(out, StencilDiagonalMatrix)
+            assert out.domain is V
+            assert out.codomain is W
+            data = out._data
+
+        # Calculate entries, or set `out=self` in default case
+        if inverse:
+            data = np.divide(1, diag, out=data)
+        elif out:
+            np.copyto(data, diag)
+        else:
+            out = self
+
+        # If needed create a new StencilDiagonalMatrix object
+        if out is None:
+            out = StencilDiagonalMatrix(V, W, data)
+
+        return out
+
 #===============================================================================
 # TODO [YG, 28.01.2021]:
 # - Check if StencilMatrix should be subclassed
diff --git a/psydac/linalg/tests/test_linalg.py b/psydac/linalg/tests/test_linalg.py
index d929e42cf..8b30802b8 100644
--- a/psydac/linalg/tests/test_linalg.py
+++ b/psydac/linalg/tests/test_linalg.py
@@ -623,22 +623,35 @@ def test_inverse_transpose_interaction(n1, n2, p1, p2, P1=False, P2=False):
 
     tol = 1e-5
     C = inverse(B, 'cg', tol=tol)
+    P = B.diagonal(inverse=True)
+
+    B_T = B.T
+    C_T = C.T
 
     # -1,T & T,-1 -> equal
-    assert isinstance(C.T, ConjugateGradient)
-    assert isinstance(inverse(B.T, 'cg', tol=tol), ConjugateGradient)
-    assert np.sqrt(sum(((C.T.dot(u) - inverse(B.T, 'cg', tol=tol).dot(u)).toarray())**2)) < 2*tol
+    assert isinstance(C_T, ConjugateGradient)
+    assert isinstance(inverse(B_T, 'cg', tol=tol), ConjugateGradient)
+    diff = C_T @ u - inverse(B_T, 'cg', tol=tol) @ u
+    assert diff.dot(diff) == 0
+
     # -1,T,T -> equal -1
-    assert np.sqrt(sum(((C.T.T.dot(u) - C.dot(u)).toarray())**2)) < 2*tol
+    diff = C_T.T @ u - C @ u
+    assert diff.dot(diff) == 0
+
     # -1,T,-1 -> equal T
-    assert isinstance(inverse(C.T, 'bicg'), BlockLinearOperator)
-    assert np.array_equal(inverse(C.T, 'bicg').dot(u).toarray(), B.T.dot(u).toarray())
+    assert isinstance(inverse(C_T, 'bicg'), BlockLinearOperator)
+    diff = inverse(C_T, 'bicg') @ u - B_T @ u
+    assert diff.dot(diff) == 0
+
     # T,-1,-1 -> equal T
-    assert isinstance(inverse(inverse(B.T, 'cg', tol=tol), 'pcg', pc='jacobi'), BlockLinearOperator)
-    assert np.array_equal(inverse(inverse(B.T, 'cg', tol=tol), 'pcg', pc='jacobi').dot(u).toarray(), B.T.dot(u).toarray())
+    assert isinstance(inverse(inverse(B_T, 'cg', tol=tol), 'pcg', pc=P), BlockLinearOperator)
+    diff = inverse(inverse(B_T, 'cg', tol=tol), 'pcg', pc=P) @ u - B_T @ u
+    assert diff.dot(diff) == 0
+
     # T,-1,T -> equal -1
-    assert isinstance(inverse(B.T, 'cg', tol=tol).T, ConjugateGradient)
-    assert np.sqrt(sum(((inverse(B.T, 'cg', tol=tol).dot(u) - C.dot(u)).toarray())**2)) < tol
+    assert isinstance(inverse(B_T, 'cg', tol=tol).T, ConjugateGradient)
+    diff = inverse(B_T, 'cg', tol=tol) @ u - C @ u
+    assert diff.dot(diff) == 0
 
     ###
     ### StencilMatrix Transpose - Inverse Tests
@@ -647,22 +660,35 @@ def test_inverse_transpose_interaction(n1, n2, p1, p2, P1=False, P2=False):
 
     tol = 1e-5
     C = inverse(S, 'cg', tol=tol)
+    P = S.diagonal(inverse=True)
+
+    S_T = S.T
+    C_T = C.T
 
     # -1,T & T,-1 -> equal
-    assert isinstance(C.T, ConjugateGradient)
-    assert isinstance(inverse(S.T, 'cg', tol=tol), ConjugateGradient)
-    assert np.sqrt(sum(((C.T.dot(v) - inverse(S.T, 'cg', tol=tol).dot(v)).toarray())**2)) < 2*tol
+    assert isinstance(C_T, ConjugateGradient)
+    assert isinstance(inverse(S_T, 'cg', tol=tol), ConjugateGradient)
+    diff = C_T @ v - inverse(S_T, 'cg', tol=tol) @ v
+    assert diff.dot(diff) == 0
+
     # -1,T,T -> equal -1
-    assert np.sqrt(sum(((C.T.T.dot(v) - C.dot(v)).toarray())**2)) < 2*tol
+    diff = C_T.T @ v - C @ v
+    assert diff.dot(diff) == 0
+
     # -1,T,-1 -> equal T
-    assert isinstance(inverse(C.T, 'bicg'), StencilMatrix)
-    assert np.array_equal(inverse(C.T, 'bicg').dot(v).toarray(), S.T.dot(v).toarray())
+    assert isinstance(inverse(C_T, 'bicg'), StencilMatrix)
+    diff = inverse(C_T, 'bicg') @ v - S_T @ v
+    assert diff.dot(diff) == 0
+
     # T,-1,-1 -> equal T
-    assert isinstance(inverse(inverse(S.T, 'cg', tol=tol), 'pcg', pc='jacobi'), StencilMatrix)
-    assert np.array_equal(inverse(inverse(S.T, 'cg', tol=tol), 'pcg', pc='jacobi').dot(v).toarray(), S.T.dot(v).toarray())
+    assert isinstance(inverse(inverse(S_T, 'cg', tol=tol), 'pcg', pc=P), StencilMatrix)
+    diff = inverse(inverse(S_T, 'cg', tol=tol), 'pcg', pc=P) @ v - S_T @ v
+    assert diff.dot(diff) == 0
+
     # T,-1,T -> equal -1
-    assert isinstance(inverse(S.T, 'cg', tol=tol).T, ConjugateGradient)
-    assert np.sqrt(sum(((inverse(S.T, 'cg', tol=tol).dot(v) - C.dot(v)).toarray())**2)) < tol
+    assert isinstance(inverse(S_T, 'cg', tol=tol).T, ConjugateGradient)
+    diff = inverse(S_T, 'cg', tol=tol) @ v - C @ v
+    assert diff.dot(diff) == 0
 
 #===============================================================================
 @pytest.mark.parametrize('n1', [3, 5])
@@ -833,8 +859,8 @@ def test_operator_evaluation(n1, n2, p1, p2):
 
     S_cg = inverse(S, 'cg', tol=tol)
     B_cg = inverse(B, 'cg', tol=tol)
-    S_pcg = inverse(S, 'pcg', pc='jacobi', tol=tol)
-    B_pcg = inverse(B, 'pcg', pc='jacobi', tol=tol)
+    S_pcg = inverse(S, 'pcg', pc=S.diagonal(inverse=True), tol=tol)
+    B_pcg = inverse(B, 'pcg', pc=B.diagonal(inverse=True), tol=tol)
     S_bicg = inverse(S, 'bicg', tol=tol)
     B_bicg = inverse(B, 'bicg', tol=tol)
     S_lsmr = inverse(S, 'lsmr', tol=tol)
@@ -944,31 +970,29 @@ def test_x0update(solver):
     # Create Inverse
     tol = 1e-6
     if solver == 'pcg':
-        A_inv = inverse(A, solver, pc='jacobi', tol=tol)
+        A_inv = inverse(A, solver, pc=A.diagonal(inverse=True), tol=tol)
     else:
         A_inv = inverse(A, solver, tol=tol)
-    options = A_inv.options
-    x0_init = options["x0"]
+
     # Check whether x0 is not None
+    x0_init = A_inv.get_options("x0")
     assert x0_init is not None
+
     # Apply inverse and check x0
     x = A_inv @ b
-    options = A_inv.options
-    x0_new1 = options["x0"]
+    x0_new1 = A_inv.get_options("x0")
     assert x0_new1 is x0_init
+
     # Change x0, apply A_inv and check for x0
     A_inv.set_options(x0 = b)
-    options = A_inv.options
-    assert options["x0"] is b
+    assert A_inv.get_options("x0") is b
+
     x = A_inv @ b
-    options = A_inv.options
-    x0_new2 = options["x0"]
-    assert x0_new2 is b
+    assert A_inv.get_options("x0") is b
+
     # Apply inverse using out=x0 and check for updated x0
-    x = A_inv.dot(b, out=x0_new2)
-    options = A_inv.options
-    x0_new3 = options["x0"]
-    assert x0_new3 is x
+    x = A_inv.dot(b, out=b)
+    assert A_inv.get_options('x0') is x
 
 #===============================================================================
 # SCRIPT FUNCTIONALITY
diff --git a/psydac/linalg/tests/test_solvers.py b/psydac/linalg/tests/test_solvers.py
index 803877257..6d4b37e73 100644
--- a/psydac/linalg/tests/test_solvers.py
+++ b/psydac/linalg/tests/test_solvers.py
@@ -59,7 +59,7 @@ def define_data(n, p, matrix_data, dtype=float):
 @pytest.mark.parametrize( 'n', [5, 10, 13] )
 @pytest.mark.parametrize('p', [2, 3])
 @pytest.mark.parametrize('dtype', [float, complex])
-@pytest.mark.parametrize('solver', ['cg', 'pcg', 'bicg', 'bicgstab', 'minres', 'lsmr', 'gmres'])
+@pytest.mark.parametrize('solver', ['cg', 'pcg', 'bicg', 'bicgstab', 'pbicgstab', 'minres', 'lsmr', 'gmres'])
 
 def test_solver_tridiagonal(n, p, dtype, solver, verbose=False):
 
@@ -67,11 +67,15 @@ def test_solver_tridiagonal(n, p, dtype, solver, verbose=False):
     # PARAMETERS
     #---------------------------------------------------------------------------
 
-    if solver in ['bicg', 'bicgstab', 'lsmr']:
+    if solver in ['bicg', 'bicgstab', 'pbicgstab', 'lsmr']:
         if dtype==complex:
             diagonals = [1-10j,6+9j,3+5j]
         else:
             diagonals = [1,6,3]
+            
+        if solver == 'pbicgstab' and dtype == complex:
+            # pbicgstab only works for real matrices
+            return
     elif solver == 'gmres':
         if dtype==complex:
             diagonals = [-7-2j,-6-2j,-1-10j]
@@ -102,10 +106,14 @@ def test_solver_tridiagonal(n, p, dtype, solver, verbose=False):
         print()
 
     #Create the solvers
-    solv  = inverse(A, solver, tol=tol, verbose=False, recycle=True)
+    if solver in ['pcg', 'pbicgstab']:
+        pc = A.diagonal(inverse=True)
+        solv = inverse(A, solver, pc=pc, tol=1e-13, verbose=verbose, recycle=True)
+    else:
+        solv = inverse(A, solver, tol=1e-13, verbose=verbose, recycle=True)
     solvt = solv.transpose()
     solvh = solv.H
-    solv2 = inverse(A@A, solver, tol=1e-13, verbose=True, recycle=True) # Test solver of composition of operators
+    solv2 = inverse(A@A, solver, tol=1e-13, verbose=verbose, recycle=True) # Test solver of composition of operators
 
     # Manufacture right-hand-side vector from exact solution
     be  = A @ xe
@@ -152,7 +160,6 @@ def test_solver_tridiagonal(n, p, dtype, solver, verbose=False):
     if solver != 'pcg':
         bc = A @ A @ xc
 
-    
     err = b - be
     err_norm = np.linalg.norm( err.toarray() )
     err2 = b2 - be2
diff --git a/psydac/linalg/utilities.py b/psydac/linalg/utilities.py
index b486bb8e6..515cef8c9 100644
--- a/psydac/linalg/utilities.py
+++ b/psydac/linalg/utilities.py
@@ -1,12 +1,18 @@
 # coding: utf-8
 
 import numpy as np
-from math                  import sqrt
+from math import sqrt
+
 from psydac.linalg.stencil import StencilVectorSpace, StencilVector
 from psydac.linalg.block   import BlockVector, BlockVectorSpace
 
-__all__ = ('array_to_psydac', 'petsc_to_psydac', '_sym_ortho')
+__all__ = (
+    'array_to_psydac',
+    'petsc_to_psydac',
+    '_sym_ortho'
+)
 
+#==============================================================================
 def array_to_psydac(x, Xh):
     """ converts a numpy array to StencilVector or BlockVector format"""
 
@@ -47,6 +53,7 @@ def array_to_psydac(x, Xh):
     u.update_ghost_regions()
     return u
 
+#==============================================================================
 def petsc_to_psydac(vec, Xh):
     """ converts a petsc Vec object to a StencilVector or a BlockVector format.
         We gather the petsc global vector in all the processes and extract the chunk owned by the Psydac Vector.
@@ -129,6 +136,7 @@ def petsc_to_psydac(vec, Xh):
     u.update_ghost_regions()
     return u
 
+#==============================================================================
 def _sym_ortho(a, b):
     """
     Stable implementation of Givens rotation.
diff --git a/psydac/utilities/utils.py b/psydac/utilities/utils.py
index 623120569..279b1bd60 100644
--- a/psydac/utilities/utils.py
+++ b/psydac/utilities/utils.py
@@ -3,8 +3,31 @@
 # Copyright 2018 Yaman Güçlü
 
 import numpy as np
+from numbers import Number
 
-__all__ = ('refine_array_1d', 'unroll_edges', 'split_space', 'split_field', 'animate_field')
+__all__ = (
+    'refine_array_1d',
+    'unroll_edges',
+    'split_space',
+    'split_field',
+    'animate_field'
+)
+
+#==============================================================================
+def is_real(x):
+    """Determine whether the given input represents a real number.
+
+    Parameters
+    ----------
+    x : Any
+
+    Returns
+    -------
+    bool
+        True if x is real, False otherwise.
+
+    """
+    return isinstance(x, Number) and np.isrealobj(x) and not isinstance(x, bool)
 
 #===============================================================================
 def refine_array_1d(x, n, remove_duplicates=True):