Inversion of a SumMatrix

NGSolve
1

I am wondering for the matrix type SumMatrix, do we have a way to take inversion?

For example, I have
hat = A.mat - B.mat@M.mat.Inverse(fesQ.FreeDofs(),inverse=“umfpack”)@B.mat.T

Do we have
“hat.Inverse”?

Thank you!

2

Which method do you want to use for the inverse? I see two black-box possibilities:

  • Use an iterative solver (CGSolver)
  • Convert the sum to a dense matrix (hat.ToDense()), and use a dense inverse

Don’t know if you like any of them. H-matrices could be nice option, but not available, currently.

3

Thanks for your response. I prefer to convert the sum to a dense matrix and use a dense inverse. I will try and update later.

4

Hi Joachim,

Is this option now available?
Trying to convert to a dense matrix crashes my notebook even for small problems.
I also tried using the CGSolver, but I don’t understand how it works.
Here is my code.

import numpy as np
from ngsolve import *
from ngsolve.webgui import Draw
from netgen.occ import *

mu0 = 4 * pi * 1e-7
eps0 = 8.85418782e-12
c0 = 299792458

rect = Rectangle(200e-3, 300e-3).Face()
rect.edges[0].name = "bottom"
rect.edges[1].name = "right"
rect.edges[2].name = "top"
rect.edges[3].name = "left"

geo = OCCGeometry(rect, dim=2)
mesh = Mesh(geo.GenerateMesh(maxh=5e-3))
# Draw(mesh)

r = y
m = 0

def a1(arg1, arg2):
    return r * curl(arg1)*curl(arg2)*dx
def a2(arg1, arg2):
    return 1/r*arg1*arg2*dx
def a3(arg1, arg2):
    return (arg1*grad(arg2) + 1/r*arg1[1]*arg2)*dx
def a4(arg1, arg2):
    return 1/r*(r**2*grad(arg1)*grad(arg2) + r*grad(arg1)[1]*arg2 + r*arg1*grad(arg2)[1] + arg1*arg2)*dx
def b1(arg1, arg2):
    return r*arg1*arg2*dx
def b2(arg1, arg2):
    return r*arg1*arg2*dx

order = 3
VrVz = HCurl(mesh, order=order, dirichlet="top")   # (E_r,E_z)
um, vm = VrVz.TnT()
G, Vphi = VrVz.CreateGradient() #Vphi -> H1
ug, vg = Vphi.TnT()
V     = VrVz * Vphi

(u, u_phi) = V.TrialFunction()
(v, v_phi) = V.TestFunction()

a = BilinearForm(a1(u, v) + a4(u_phi, v_phi) + m**2*a2(u, v) + m*a3(u, v_phi) + m*a3(v, u_phi))
b = BilinearForm(b1(u, v) + b2(u_phi, v_phi))

mc = BilinearForm(b1(um, vm))
mg = BilinearForm(b2(ug, vg))

apre = BilinearForm(a1(u, v) + a4(u_phi, v_phi) + m**2*a2(u, v) + m*a3(u, v_phi) + m*a3(v, u_phi) + b1(u, v) + b2(u_phi, v_phi))
pre = Preconditioner(apre, type="direct", inverse="sparsecholesky")

with TaskManager():
    a.Assemble()
    b.Assemble()
    mc.Assemble()
    mg.Assemble()
    apre.Assemble()
    
    I_VrVz = IdentityMatrix(VrVz.ndof)
    I_phi = IdentityMatrix(Vphi.ndof)
    
    GT = G.CreateTranspose()
    S = GT @ mc.mat @ G + mg.mat
    Sinv = S.Inverse(inverse="sparsecholesky", freedofs=Vphi.FreeDofs()) #<--Inverse of SumMatrix, does not work
#     Sinv = solvers.CGSolver(S, I_phi) <--- Tried CGSolver, not sure how to build a preconditioner for S

    embHc, embH1 = V.embeddings
    proj =  embHc @  (I_VrVz - G@Sinv@GT@mc.mat) @ embHc.T +\
    embH1 @ (I_phi - Sinv@mg.mat) @ embH1.T +\
    embH1 @ (-Sinv@GT@mc.mat) @ embHc.T +\
    embHc @ (-G@Sinv@mg.mat) @ embH1.T
#     print (proj.GetOperatorInfo())
    projpre = proj @ pre.mat

    evals, evecs = solvers.PINVIT(a.mat, b.mat, pre=projpre, num=5, maxit=20,
                                  printrates=True)

Best regards