From 578bf589f7d99e5a08683f2515951edab5f9b0c7 Mon Sep 17 00:00:00 2001 From: Umberto Zerbinati Date: Wed, 19 Jun 2024 10:22:11 +0100 Subject: [PATCH] Solver now can be an operator in matrix free form Signed-off-by: Umberto Zerbinati --- ngsPETSc/ksp.py | 45 +++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 43 insertions(+), 2 deletions(-) diff --git a/ngsPETSc/ksp.py b/ngsPETSc/ksp.py index 411a754..6f1e4e6 100644 --- a/ngsPETSc/ksp.py +++ b/ngsPETSc/ksp.py @@ -4,7 +4,7 @@ ''' from petsc4py import PETSc -from ngsolve import la, GridFunction +from ngsolve import la, GridFunction, BaseMatrix from ngsPETSc import Matrix, VectorMapping @@ -45,7 +45,7 @@ def __init__(self, a, fes, p=None, solverParameters=None, optionsPrefix=None, nu for optName, optValue in solverParameters.items(): options_object[optName] = optValue - #Creating the PETSc Matrix + #Creating the PETSc Matrix A = Matrix(Amat, fes).mat A.setOptionsPrefix(optionsPrefix) A.setFromOptions() @@ -90,3 +90,44 @@ def view(self): ''' self.ksp.view() + + def asOperator(self, ngsMat): + return self.Operator(self.ksp, ngsMat) + + class Operator(BaseMatrix): + def __init__(self, ksp, ngsMat): + self.ksp = ksp + self.ngsMat = ngsMat + def Shape(self): + ''' + Shape of the BaseMatrix + + ''' + return self.ngsMat.shape + + def CreateVector(self,col): + ''' + Create vector corresponding to the matrix + + :arg col: True if one want a column vector + + ''' + return self.ngsMat.CreateVector(not col) + + def Mult(self,x,y): + ''' + BaseMatrix multiplication Ax = y + :arg x: vector we are multiplying + :arg y: vector we are storeing the result in + + ''' + self.ksp.solve(x,y) + + def MultTrans(self,x,y): + ''' + BaseMatrix multiplication A^T x = y + :arg x: vector we are multiplying + :arg y: vector we are storeing the result in + + ''' + raise NotImplementedError("Transpose multiplication not implemented")