I would like to solve a dense linear system the form in python
$$
L\left(\boldsymbol{x}\right):=\left[\gamma^+\left[\boldsymbol{A}+\frac{1}{2}\boldsymbol{B}^{-1}\right]
+\gamma^-\left[\boldsymbol{A}-\frac{1}{2}\boldsymbol{B}^{-1}\right]\right]\cdot\boldsymbol{x}=\boldsymbol{b}
$$
I thought it is a good idea not to explicitly compute the inverse of $\boldsymbol{B}$. Therefore, I implemented the Operator $L$ to solve $L\left(\boldsymbol{x}\right)=\boldsymbol{b}$ using Krylov methods like cg or gmres. I am using the scipy.sparse.linalg.LinearOperator
class for the operator see the docs here. The product $\boldsymbol{B}^{-1}\cdot\boldsymbol{x}$ computed by another Krylov iteration or an L-U decomposition depending on system size.
However, for larger problems I would like to improve the rate of convergence of the outer iteration. I neither have sparse matrices nor an explicit representation of my matrix. Therefore, as far as I know the classical preconditioners for Krylov methods like ilu or Jacobi's method are not applicable.
Are there other methods which can be used? And are there python libraries for these methods?