I have tried implementing the MHSS preconditioner for the CG solver in NGSolve.

P = (1 + j)(B + C) where B is the real part and C is the imaginary part.
I could calculate a_pre_inv = a_pre.Inverse(fes_pre.FreeDofs(), inverse="sparsecholesky") where a_pre is B + C.
Is there any way to scale this inverse with 1/(1-1j) and use this matrix as a preconditioner for the CG solver of the complex system?

Thank you for your response. I have time to look into it again.
However, the solver did not converge. I have implemented this preconditioner in Scipy by exporting the matrix from NGSolve and it does work. When exporting the matrix, I removed the row and column with fixed DOFs. I reckon that the CG solver did not converge because the fixed DOFs are also taken into consideration if I just used:

I figured it out. when constructing the bilinear form of a_pre, I set the complex = False for H1 because it was a real number space. I tried with complex = True. It works. The solver does converge.

NGSolve contains a built-in scaled matrix, you can write s*matrix to generate a scaled matrix.
And, there is a real2complex wrapper, which allows to apply a real matrix to complex vectors (by applying it to the real, and the imaginary parts separately):

from ngsolve.la import Real2ComplexMatrix
Ac = 1/(1-1j) * Real2ComplexMatrix(A)

Thank you for your suggestion. Indeed, the inverse is now faster as it is a real space. I changed the code to

class ScaledBaseMat(BaseMatrix):
def __init__(self, mat):
super().__init__()
self.mat = mat
def Mult(self, x: BaseVector, y: BaseVector):
Ac = 1/(1+1j) * Real2ComplexMatrix(self.mat)
y.data = Ac * x

It works for direct inverse: a_pre_inv = a_pre_b.mat.Inverse(fes_pre.FreeDofs() , inverse="sparsecholesky")

When dealing with a bigger matrix, this inverse is expensive, therefore, I used iterative solver to compute:

can you send a MWE for that ? It looks like we miss to keep a shared object somewhere.

You can directly use Ac = 1/(1+1j) * Real2ComplexMatrix(self.mat) instead of your own ScaledBaseMat. Justs create it once, and keep that BaseMatrix (= linear operator). Does it lead to the same problem ?

I implemented your suggestion, but it has not worked.
I think about the problem again. Even if one defines a_pre_inv = CGSolver(a_pre_b.mat, pre=preJpoint) as a real matrix and then computed y.data = 1/(1+1j) * a_pre_inv * x (applying a preconditioner to a vector x). The complex problem is solved as x is a complex vector.

This can be organized as the solution of two real systems. One for the real part and one for the imaginary part. I have tried:

However, it does not work. I have successfully implemented the BiCGStab solver in NGSolve. Any suggestion to improve the BiCGStab solver and apply the MHSS preconditioner using an iterative solver is highly appreciated.

I attached the MWE script. Thank you so much for your time.

Is there any way to set a real vector to a real part of a complex vector and another real vector to an imaginary part of the same complex vector? The workaround below gives the expected results, but there might be a more efficient way to avoid multiplication.