Hello,
I need to solve a non-symmetric linear system and then later also solve the adjoint system. I tried to achieve this by just factorizing the system matrix once using umfpack and then using the transpose via the .T option. Unfortunately, this does not seem to work correctly as this simplified example shows:
from ngsolve import *
from ngsolve.meshes import Make1DMesh
from ngsolve.la import BaseVector
import numpy as np
mesh=Make1DMesh(n=10)
fes = H1(mesh, order=1, dirichlet=".*")
N=fes.ndof
u,v=fes.TnT()
bf = BilinearForm(fes)
bf += grad(u)*grad(v)*dx-u*grad(v)*dx#non-symmetric bilinear form
bf.Assemble()
print(bf.mat.GetInverseType())#umfpack
Binv=bf.mat.Inverse(fes.FreeDofs())
a=BaseVector(np.linspace(0,10,N))
b=BaseVector(np.linspace(-1,1,N))
Binv_a=BaseVector(N)
BinvT_b=BaseVector(N)
BTinv_b=BaseVector(N)
Binv_a.data=Binv*a
BinvT_b.data=Binv.T*b #Does not work correctly
BTinv_b.data=bf.mat.CreateTranspose().Inverse(fes.FreeDofs())*b #Does work correctly
print(InnerProduct(Binv_a,b))#3.8729422763480676
print(InnerProduct(a,BinvT_b))#1.2909807587826898
print(InnerProduct(a,BTinv_b))#3.8729422763480676
Is there another way to do this in ngsolve using only a single factorization?
Thank you in advance for your help.