How to use the petsc4py.PETSc.KSP function in petsc4py

A.setPythonContext(pde)
A.setUp()

y = da.createGlobalVec()

pde.formRHS(x, val=1.0)
A.mult(x, b)
pde.formRHS(y, val=-2.0 * (2.0 * np.pi) ** 2)

# u = da.createNaturalVec()
# da.globalToNatural(b, u)

# print((b - y).norm(PETSc.NormType.NORM_INFINITY))
# exit()

ksp = PETSc.KSP().create()
ksp.setType('cg')
pc = ksp.getPC()
pc.setType('ilu')
ksp.setFromOptions()

ksp.setOperators(A)

x1 = da.createGlobalVec()
pde.formRHS(b,val=1)
ksp.solve(b, x1)

x2 = da.createGlobalVec()
pde.formRHS(b,val=0)
ksp.solve(b, x2)

print(x1.array)

petsc4py.init(sys.argv)
    from petsc4py import PETSc
    OptDB = PETSc.Options()
    for k,v in precon_opts[precond_type].items():
        OptDB[k]=v
    A_petsc = PETSc.Mat().createAIJ(A.shape,
                                    csr=(A.indptr,
                                         A.indices,
                                         A.data))
    # obtain vectors for storing
    # the solution  and the rhs
    x_petsc, b_petsc = A_petsc.getVecs()
    # fill the rhs PETSc vector
    # from the rhs numpy array
    b_petsc[...] = RHS # note the syntax sugar
    ksp = PETSc.KSP().create()
    ksp.setFromOptions()  # configure from OptDB
    
    ksp.setInitialGuessNonzero(guess)
    ksp.setOperators(A_petsc)   
    ksp.setTolerances(1e-10)
    x_petsc[...] = cur_val
    ksp.solve(b_petsc, x_petsc) 
    x_petsc[...] = cur_val
    gc.collect()
    t.sleep(1)
    t1 = t.time()
    ksp.solve(b_petsc, x_petsc) 
    t2 = t.time()
    x_petsc[...] = cur_val
    t3 = t.time()
    ksp.solve(b_petsc, x_petsc)

# invoke super init, passing number of dofs, dtype_u and dtype_f
        super(petsc_grayscott_multiimplicit, self).__init__(init=da, dtype_u=dtype_u, dtype_f=dtype_f,
                                                            params=problem_params)

        # compute dx, dy and get local ranges
        self.dx = 100.0 / (self.params.nvars[0])
        self.dy = 100.0 / (self.params.nvars[1])
        (self.xs, self.xe), (self.ys, self.ye) = self.init.getRanges()

        # compute discretization matrix A and identity
        self.A = self.__get_A()
        self.Id = self.__get_Id()
        self.localX = self.init.createLocalVec()

        # setup linear solver
        self.ksp = PETSc.KSP()
        self.ksp.create(comm=self.params.comm)
        self.ksp.setType('cg')
        pc = self.ksp.getPC()
        pc.setType('none')
        self.ksp.setInitialGuessNonzero(True)
        self.ksp.setFromOptions()
        self.ksp.setTolerances(rtol=self.params.lsol_tol, atol=self.params.lsol_tol, max_it=self.params.lsol_maxiter)
        self.ksp_itercount = 0
        self.ksp_ncalls = 0

        # setup nonlinear solver
        self.snes = PETSc.SNES()
        self.snes.create(comm=self.params.comm)
        # self.snes.getKSP().setType('cg')
        # self.snes.setType('ngmres')
        self.snes.setFromOptions()

if self.have_MPI == True:
            gamma_array_DispX = PETSc.Vec().create(self.comm)
            gamma_array_DispY = PETSc.Vec().create(self.comm)
            gamma_array_DispZ = PETSc.Vec().create(self.comm)
            gamma_array_DispX.setType('mpi')
            gamma_array_DispY.setType('mpi')
            gamma_array_DispZ.setType('mpi')
            KSP_solver = PETSc.KSP().create(self.comm)
          else:
            gamma_array_DispX = PETSc.Vec().create()
            gamma_array_DispY = PETSc.Vec().create()
            gamma_array_DispZ = PETSc.Vec().create()
            gamma_array_DispX.setType('seq')
            gamma_array_DispY.setType('seq')
            gamma_array_DispZ.setType('seq')
            KSP_solver = PETSc.KSP().create()
          gamma_array_DispX.setSizes(self.nSolidInterfacePhysicalNodes+self.d_RBF)
          gamma_array_DispY.setSizes(self.nSolidInterfacePhysicalNodes+self.d_RBF)
          gamma_array_DispZ.setSizes(self.nSolidInterfacePhysicalNodes+self.d_RBF)
          gamma_array_DispX.set(0.0)
          gamma_array_DispY.set(0.0)
          gamma_array_DispZ.set(0.0)
          KSP_solver.setType('fgmres')
          KSP_solver.getPC().setType('jacobi')
          KSP_solver.setOperators(self.MappingMatrixA)
          KSP_solver.setFromOptions()
          #print(KSP_solver.getInitialGuessNonzero())
          KSP_solver.setInitialGuessNonzero(True)
          #print(KSP_solver.getInitialGuessNonzero())
          KSP_solver.solve(self.solidInterface_array_DispX, gamma_array_DispX)
          KSP_solver.solve(self.solidInterface_array_DispY, gamma_array_DispY)
          KSP_solver.solve(self.solidInterface_array_DispZ, gamma_array_DispZ)

self.ssarun = ssarun
        self.ssa = self.ssarun.ssa
        self.grid = ssarun.grid

        self.tmpV = PISM.IceModelVec2V()
        self.tmpV.create(self.grid, "work vector (2V)", PISM.WITHOUT_GHOSTS, WIDE_STENCIL)

        self.tmpS = PISM.IceModelVec2S()
        self.tmpS.create(self.grid, "work vector (2S)", PISM.WITHOUT_GHOSTS, WIDE_STENCIL)

        self.tmpS2 = PISM.IceModelVec2S()
        self.tmpS2.create(self.grid, "work vector (2S)", PISM.WITHOUT_GHOSTS, WIDE_STENCIL)

        ksp_rtol = 1e-12
        self.ksp = PETSc.KSP()
        self.ksp.create(self.grid.com)
        self.ksp.setTolerances(ksp_rtol, PETSc.DEFAULT, PETSc.DEFAULT, PETSc.DEFAULT)
        self.ksp.getPC().setType('bjacobi')
        self.ksp.setFromOptions()

        (self.designFunctional, self.stateFunctional) = PISM.invert.ssa.createGradientFunctionals(ssarun)

        self.designForm = None

def _solve(self,X,rhs,Npts,Nctl):
        '''Solve for the control points'''
        if not self.NO_PRINT:
            print 'LMS solving...'

        if USE_PETSC:
            ksp = PETSc.KSP()
            ksp.create(PETSc.COMM_WORLD)
            ksp.getPC().setType('none')
            ksp.setType('lsqr')
            #ksp.setInitialGuess(True)
            print 'Iteration   Residual'
            def monitor(ksp, its, rnorm):
                if mod(its,100) == 0:
                    print '%5d      %20.15g'%(its,rnorm)

            ksp.setMonitor(monitor)
            ksp.setTolerances(rtol=1e-15, atol=1e-15, divtol=100, max_it=500)

            ksp.setOperators(self.J)
            ksp.solve(rhs, X) 

            coef_temp = zeros((Nctl,3))

da = da.refine()

        # invoke super init, passing number of dofs, dtype_u and dtype_f
        super(heat2d_petsc_forced, self).__init__(init=da, dtype_u=dtype_u, dtype_f=dtype_f, params=problem_params)

        # compute dx, dy and get local ranges
        self.dx = 1.0 / (self.init.getSizes()[0] - 1)
        self.dy = 1.0 / (self.init.getSizes()[1] - 1)
        (self.xs, self.xe), (self.ys, self.ye) = self.init.getRanges()

        # compute discretization matrix A and identity
        self.A = self.__get_A()
        self.Id = self.__get_Id()

        # setup solver
        self.ksp = PETSc.KSP()
        self.ksp.create(comm=self.params.comm)
        self.ksp.setType('gmres')
        pc = self.ksp.getPC()
        pc.setType('none')
        # self.ksp.setInitialGuessNonzero(True)
        self.ksp.setFromOptions()
        self.ksp.setTolerances(rtol=self.params.sol_tol, atol=self.params.sol_tol, max_it=self.params.sol_maxiter)

        self.ksp_ncalls = 0
        self.ksp_itercount = 0

# Want the iterations to start from random (iterative)
    wh.block_vec().randomize()

    # Default is minres
    if '-ksp_type' not in solver_params: solver_params['-ksp_type'] = 'minres'
            
    opts = PETSc.Options()
    for key, value in solver_params.iteritems():
        opts.setValue(key, None if value == 'none' else value)
        
    ksp = PETSc.KSP().create()
    ksp.setOperators(ii_PETScOperator(AA))
    ksp.setPC(ii_PETScPreconditioner(BB, ksp))
    
    ksp.setNormType(PETSc.KSP.NormType.NORM_PRECONDITIONED)
    ksp.setFromOptions()
    
    ksp.solve(as_petsc_nest(bb), wh.petsc_vec())
    
    niters = ksp.getIterationNumber()

    return wh, niters

assert isinstance(par_L,ParMat_petsc4py)
        
        self.pccontext = None
        self.preconditioner = None
        self.pc = None
        self.solverName  = "PETSc"
        self.par_fullOverlap = True
        self.par_firstAssembly=True
        self.par_L   = par_L
        self.petsc_L = par_L
        self.csr_rep_local = self.petsc_L.csr_rep_local
        self.csr_rep = self.petsc_L.csr_rep
        self.bdyNullSpace = bdyNullSpace

        # create petsc4py KSP object and attach operators
        self.ksp = p4pyPETSc.KSP().create()
        self._setMatOperators()
        self.ksp.setOperators(self.petsc_L,self.petsc_L)

        self.setResTol(rtol_r,atol_r)

        convergenceTest = 'r-true'
        if convergenceTest == 'r-true':
            self.r_work = self.petsc_L.getVecLeft()
            self.rnorm0 = None
            self.ksp.setConvergenceTest(self._converged_trueRes)
        else:
            self.r_work = None

        if prefix is not None:
            self.ksp.setOptionsPrefix(prefix)
        if Preconditioner is not None: