Symmetrized Gauss-Seidel preconditioner in Parallel (MPI)

General
1

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

2

no, there is no MPI-parallel GS in NGSolve.
We recommend to use ngspetsc for that.

Joachim

3

Thanks Joachim for the clarification.

Best
Philipp

4

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