Hello,

I am fairly new to preconditioning large algebraic systems and facing problems trying to construct a symmetrized Gauss-Seidel preconditioner, that works in parallel with MPI.

The following Code works without MPI.

```
## Symmetrized GS preconditioner ###
class SymmetricGS(BaseMatrix):
def __init__ (self, smoother):
super(SymmetricGS, self).__init__()
self.smoother = smoother
def Mult (self, x, y):
y[:] = 0.0
self.smoother.Smooth(y, x)
self.smoother.SmoothBack(y,x)
def Height (self):
return self.smoother.height
def Width (self):
return self.smoother.height
preJpoint = a.mat.CreateSmoother(fes.FreeDofs(True))
preGS = SymmetricGS(preJpoint)
```

Is there a way to construct a symmetrized GS preconditioner, that works with MPI?

Best

Philipp

no, there is no MPI-parallel GS in NGSolve.

We recommend to use ngspetsc for that.

Joachim

Thanks Joachim for the clarification.

Best

Philipp

Hello Joachim and others,

unfortunately there seems to be quite a major Bug within ngsPETSc, when i want to wrap PETSc preconditioner, and wrap it into an NGSolve preconditioner. This holds for nonlinear problems, where the preconditioner needs to be updated within a newton raphson scheme.

When i use the wrapper from ngsPETSc this update is not done and the algorithm simply stops after the first iteration.

Here is a simple Newton solver which showcases the problem:

```
### Simple Newton Solver ###
from ngsolve.solvers import *
from ngsolve.la import EigenValues_Preconditioner
from ngsPETSc.pc import *
def SimpleNewtonSolve(a,gfsol,tol=1e-6,maxits=8):
res = gfsol.vec.CreateVector()
dsol = gfsol.vec.CreateVector()
## c = Preconditioner(a,'local')
## buit in Jacobi preconditioner works fine ##
c = Preconditioner(a, "PETScPC", pc_type="jacobi")
## Jacobi Precondtioner @ PETSc is not updated
for it in range(maxits):
print ("Iteration {:3} ".format(it),end="")
a.Apply(gfsol.vec, res)
a.AssembleLinearization(gfsol.vec)
lams = EigenValues_Preconditioner(mat=a.mat, pre=c.mat)
print("condition number:", max(lams)/min(lams))
# solve linear system of equations
solvers.GMRes(A = a.mat, pre = c.mat, b = res,
x = dsol.data, printrates="\r", tol = 1e-10, maxsteps = 500)
gfsol.vec.data -= dsol
#stopping criteria
stopcritval = sqrt(abs(InnerProduct(dsol,res)))
print ("<A u",it,", A u",it,">_{-1}^0.5 = ", stopcritval)
if stopcritval < tol:
break
```

Best

Philipp