Scale the inverse of matrix

NGSolve
1

Dear all,

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?

Any help is highly appreciated.

Best regards,
VC

2

Maybe there is even an easier way, but one way is to derive from BaseMatrix and implement the scaling there. Something like this:

class ScaledBaseMat(BaseMatrix):
  def __init__(self, mat):
    self.mat = mat
    super().__init__()

  def Mult(self, x, y):
    y.data = self.mat * x
    y[:] *= 1/(1-1j)

and then use this matrix in your preconditioner.

See here for an example of such a BaseMatrix overload:
https://docu.ngsolve.org/latest/i-tutorials/unit-2.1.2-blockjacobi/blockjacobi.html#Implement-a-forward-backward-GS-preconditioner

1 Like
3

Dear Christopher,

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:

a_pre_inv = a_pre.Inverse(fes.FreeDofs())
pre = ScaledBaseMat(a_pre_inv)
inv = CGSolver(a.mat, pre, maxiter=2000, tol=1e-8, printrates=True, conjugate = False)

I wondered if there should be some way to deal with fixed DOFs when creating a Preconditioner in NGSolve.
Do you have any suggestions?

Thank you so much for your time.

Best regards,
VC

4

Dear Christopher,

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.

Best regards,
VC

5

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)

best, Joachim

6

Dear Jachim,

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:

preJpoint = a_pre_b.mat.CreateSmoother(fes_pre.FreeDofs())
a_pre_inv = CGSolver(a_pre_b.mat, pre=preJpoint, maxiter=1000, tol=1e-2, printrates=False)

However, I got the error when executing y.data = Ac * x
RuntimeError: bad_weak_ptr

do you have any suggestions about it?

Best regards,
VC

8

Hi VienC,

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 ?

Joachim

9

Dear Joachim,

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:

class ScaledBaseMat(BaseMatrix):
  def __init__(self, mat, freedofs: BitArray = None):      
    super().__init__()
    self.mat = mat

  def Mult(self, x, y):

    x_real = x.FV().NumPy().real
    x_imag = x.FV().NumPy().imag

    x1 = BaseVector(0.5*(x_real + x_imag))
    x2 = BaseVector(0.5*(x_imag - x_real))
    y_real = self.mat * x1
    y_imag = self.mat * x2
    y.data =  y_real + 1j * y_imag

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.

Best regards,
VC
MHSS_iterative.py (14.6 KB)

10

Dear Joachim,

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.

y.data = -1j*1j*y_real
y.data += 1j*y_imag

Best regards,
Vien Che

11

with the .FV() you get a view of the vector memory, with it you can extract real/imaginary components.

v = CreateVVector(5, complex=True)
v[:] = 1+2j

vr = CreateVVector(5)
vi = CreateVVector(5)

vr.FV()[:] = v.FV().real
vi.FV()[:] = v.FV().imag

print (v, vr, vi)