# Authors: # Loic Gouarin # Nicole Spillane # # License: BSD 3 clause from petsc4py import PETSc from slepc4py import SLEPc import mpi4py.MPI as mpi import numpy as np import os from .projection import projection, GenEO_V0, minimal_V0, coarse_operators class PCAWG: def __init__(self, A): print('PCAWG !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') OptDB = PETSc.Options() self.switchtoASM = OptDB.getBool('PCNew_switchtoASM', False) #use Additive Schwarz as a preconditioner instead of BNN self.switchtoASMpos = OptDB.getBool('PCNew_switchtoASMpos', False) #use Additive Schwarz as a preconditioner instead of BNN self.verbose = OptDB.getBool('PCNew_verbose', False) self.GenEO = OptDB.getBool('PCNew_GenEO', True) #self.H2addCS = OptDB.getBool('PCNew_H2addCoarseSolve', True) self.H2projCS = OptDB.getBool('PCNew_H2CoarseProjection', True) self.H3addCS = OptDB.getBool('PCNew_H3addCoarseSolve', True) self.H3projCS = OptDB.getBool('PCNew_H3CoarseProjection', True) self.compute_ritz_apos = OptDB.getBool('PCNew_ComputeRitzApos', False) self.nev = OptDB.getInt('PCNew_Bs_nev', 20) #number of vectors asked to SLEPc for cmputing negative part of Bs self.viewPC = OptDB.getBool('PCNew_view', True) self.viewV0 = OptDB.getBool('PCNew_viewV0', False) self.viewGenEOV0 = OptDB.getBool('PCNew_viewGenEO', False) self.viewminV0 = OptDB.getBool('PCNew_viewminV0', False) self.viewnegV0 = OptDB.getBool('PCNew_viewnegV0', False) self.test_case = OptDB.getString('test_case', 'default') self.H2addCS = True #OptDB.getBool('PCNew_H2addCoarseSolve', True) (it is currently not an option to use a projected preconditioner for H2) ksp = PETSc.KSP().create() ksp.setOperators(A) ksp.setOptionsPrefix("H_") pc = ksp.pc pc.setType('asm') pc.setASMOverlap(2) pc.setASMType(PETSc.PC.ASMType.BASIC) #0: none (block Jacobi); 1:restrict (non sym); 2: interpolate (non sym); 3: basic pc.setFromOptions() pc.setUp() isb, isl = pc.getASMLocalSubdomains() lgmap = PETSc.LGMap().create(isb[0]) localksp = pc.getASMSubKSP()[0] As = localksp.getOperators()[0] global_sizes = A.getSizes() local_sizes = As.getSizes() works_1, works_2 = As.getVecs() work_1, work_2 = A.getVecs() scatter_l2g = PETSc.Scatter().create(works_1, None, work_1, isb[0]) Bs = As.duplicate() for i in range(global_sizes[0][1]): work_1.set(0.) work_1[i] = 1. work_1.assemble() scatter_l2g(work_1, works_1, PETSc.InsertMode.INSERT_VALUES, PETSc.ScatterMode.SCATTER_REVERSE) As.mult(works_1, works_2) work_2.set(0.) scatter_l2g(works_2, work_2, PETSc.InsertMode.ADD_VALUES) col = np.asarray(np.where(works_1[:] == 1)[0], dtype='int32') scatter_l2g(work_2, works_1, PETSc.InsertMode.INSERT_VALUES, PETSc.ScatterMode.SCATTER_REVERSE) if col.size > 0: mask_B = As[:, col[0]] != 0 mask_works = works_1[:] != 0 mask = np.logical_and(mask_B, mask_works) indices = np.asarray(range(local_sizes[0][0]), dtype='int32') Bs.setValues(indices[mask], col[0], As[:, col[0]][mask]**2/works_1[:][mask], PETSc.InsertMode.INSERT_VALUES) Bs.assemble() #mumps solver for Bs Bs_ksp = PETSc.KSP().create(comm=PETSc.COMM_SELF) Bs_ksp.setOptionsPrefix("Bs_ksp_") Bs_ksp.setOperators(Bs) Bs_ksp.setType('preonly') Bs_pc = Bs_ksp.getPC() Bs_pc.setType('cholesky') Bs_pc.setFactorSolverType('mumps') Bs_pc.setFactorSetUpSolverType() Bs_pc.setUp() Bs_ksp.setFromOptions() mus = works_1.duplicate() #compute the multiplicity of each dof works_1.set(1.) work_1.set(0.) scatter_l2g(works_1, work_1, PETSc.InsertMode.ADD_VALUES) scatter_l2g(work_1, mus, PETSc.InsertMode.INSERT_VALUES, PETSc.ScatterMode.SCATTER_REVERSE) _, mult_max = work_1.max() if self.viewPC: _, self.ns = mus.getSizes() _, self.nglob = work_1.getSizes() tempglobal = work_1.getArray(readonly=True) templocal = mus.getArray(readonly=True) self.nints = np.count_nonzero(tempglobal == 1) #interor dofs in this subdomain self.nGammas = np.count_nonzero(templocal -1) #interor dofs in this subdomain invmus = mus.duplicate() invmus = 1/mus if mpi.COMM_WORLD.rank == 0: print(f'multmax: {mult_max}') DVnegs = [] Vnegs = [] invmusVnegs = [] #BEGIN diagonalize Bs #Eigenvalue Problem for smallest eigenvalues eps = SLEPc.EPS().create(comm=PETSc.COMM_SELF) eps.setDimensions(nev=self.nev) eps.setProblemType(SLEPc.EPS.ProblemType.HEP) eps.setOperators(Bs) #print(f'dimension of Bs : {Bs.getSize()}') #OPTION 1: works but dense algebra eps.setType(SLEPc.EPS.Type.LAPACK) eps.setWhichEigenpairs(SLEPc.EPS.Which.SMALLEST_REAL) #with lapack this just tells slepc how to order the eigenpairs ##END OPTION 1 ##OPTION 2: default solver (Krylov Schur) but error with getInertia - is there a MUMPS mattype - Need to use MatCholeskyFactor #if Which eigenpairs is set to SMALLEST_REAL, some are computed but not all ##Bs.setOption(PETSc.Mat.Option.SYMMETRIC, True) ##Bs.convert('sbaij') ##IScholBs = is_A.duplicate() ##Bs.factorCholesky(IScholBs) #not implemented #tempksp = PETSc.KSP().create(comm=PETSc.COMM_SELF) #tempksp.setOperators(Bs) #tempksp.setType('preonly') #temppc = tempksp.getPC() #temppc.setType('cholesky') #temppc.setFactorSolverType('mumps') #temppc.setFactorSetUpSolverType() #tempF = temppc.getFactorMatrix() #tempF.setMumpsIcntl(13, 1) #needed to compute intertia according to slepcdoc, inertia computation still doesn't work though #temppc.setUp() ##eps.setOperators(tempF) #eps.setWhichEigenpairs(SLEPc.EPS.Which.ALL) #eps.setInterval(PETSc.NINFINITY,0.0) #eps.setUp() ##eps.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_REAL) ##eps.setTarget(0.) ##if len(Vnegs) > 0 : ## eps.setDeflationSpace(Vnegs) ##if mpi.COMM_WORLD.rank == 0: ## eps.view() ##END OPTION 2 eps.solve() if eps.getConverged() < self.nev: PETSc.Sys.Print('for Bs in subdomain {}: {} eigenvalues converged (less that the {} requested)'.format(mpi.COMM_WORLD.rank, eps.getConverged(), self.nev), comm=PETSc.COMM_SELF) Dnegs = [] Dposs = [] for i in range(eps.getConverged()): tempscalar = np.real(eps.getEigenvalue(i)) if tempscalar < 0. : Dnegs.append(-1.*tempscalar) Vnegs.append(works_1.duplicate()) eps.getEigenvector(i,Vnegs[-1]) DVnegs.append(Dnegs[-1] * Vnegs[-1]) invmusVnegs.append(invmus * Vnegs[-1]) else : Dposs.append(tempscalar) if self.verbose: PETSc.Sys.Print('for Bs in subdomain {}: ncv= {} with {} negative eigs'.format(mpi.COMM_WORLD.rank, eps.getConverged(), len(Vnegs), self.nev), comm=PETSc.COMM_SELF) print(f'values of Dnegs {np.array(Dnegs)}') nnegs = len(Dnegs) #print(f'length of Dnegs {nnegs}') #END diagonalize Bs if self.viewnegV0: print('###') print(f'view of Vneg in Subdomain {mpi.COMM_WORLD.rank}') print(f'ncv = {eps.getConverged()} eigenvalues converged') print(f'{nnegs=}') print(f'values of Dnegs: {np.array(Dnegs)}') works_1.set(0.) RsVnegs = [] Vneg = [] Dneg = [] RsDVnegs = [] RsDnegs = [] for i in range(mpi.COMM_WORLD.size): nnegi = len(Vnegs) if i == mpi.COMM_WORLD.rank else None nnegi = mpi.COMM_WORLD.bcast(nnegi, root=i) for j in range(nnegi): Vneg.append(Vnegs[j].copy() if i == mpi.COMM_WORLD.rank else works_1.copy()) dnegi = Dnegs[j] if i == mpi.COMM_WORLD.rank else None dnegi = mpi.COMM_WORLD.bcast(dnegi, root=i) Dneg.append(dnegi) #print(f'i Dneg[i] = {i} {Dneg[i]}') for i, vec in enumerate(Vneg): work_1.set(0) scatter_l2g(vec, work_1, PETSc.InsertMode.ADD_VALUES) scatter_l2g(work_1, works_1, PETSc.InsertMode.INSERT_VALUES, PETSc.ScatterMode.SCATTER_REVERSE) if works_1.norm() != 0: RsVnegs.append(works_1.copy()) RsDVnegs.append(Dneg[i]*works_1.copy()) RsDnegs.append(Dneg[i]) #TO DO: here implement RsVnegs and RsDVnegs #self.Vneg = Vneg # self.Vnegs = Vnegs # self.DVnegs = DVnegs # self.scatterl #Local Apos and Aneg Aneg = PETSc.Mat().createPython([work_1.getSizes(), work_1.getSizes()], comm=PETSc.COMM_WORLD) Aneg.setPythonContext(Aneg_ctx(Vnegs, DVnegs, scatter_l2g, works_1, works_2)) Aneg.setUp() Apos = PETSc.Mat().createPython([work_1.getSizes(), work_1.getSizes()], comm=PETSc.COMM_WORLD) Apos.setPythonContext(Apos_ctx(A, Aneg )) Apos.setUp() #A pos = A + Aneg so it could be a composite matrix rather than Python type Anegs = PETSc.Mat().createPython([works_1.getSizes(), works_1.getSizes()], comm=PETSc.COMM_SELF) Anegs.setPythonContext(Anegs_ctx(Vnegs, DVnegs)) Anegs.setUp() Aposs = PETSc.Mat().createPython([works_1.getSizes(), works_1.getSizes()], comm=PETSc.COMM_SELF) Aposs.setPythonContext(Aposs_ctx(Bs, Anegs )) Aposs.setUp() projVnegs = PETSc.Mat().createPython([works_1.getSizes(), works_1.getSizes()], comm=PETSc.COMM_SELF) projVnegs.setPythonContext(projVnegs_ctx(Vnegs)) projVnegs.setUp() projVposs = PETSc.Mat().createPython([works_1.getSizes(), works_1.getSizes()], comm=PETSc.COMM_SELF) projVposs.setPythonContext(projVposs_ctx(projVnegs)) projVposs.setUp() #TODO Implement RsAposRsts, this is the restriction of Apos to the dofs in this subdomain. So it applies to local vectors but has non local operations RsAposRsts = PETSc.Mat().createPython([works_1.getSizes(), works_1.getSizes()], comm=PETSc.COMM_SELF) #or COMM_WORLD ? RsAposRsts.setPythonContext(RsAposRsts_ctx(As,RsVnegs,RsDVnegs)) RsAposRsts.setUp() invAposs = PETSc.Mat().createPython([works_1.getSizes(), works_1.getSizes()], comm=PETSc.COMM_SELF) invAposs.setPythonContext(invAposs_ctx(Bs_ksp, projVposs )) invAposs.setUp() ksp_Aposs = PETSc.KSP().create(comm=PETSc.COMM_SELF) ksp_Aposs.setOperators(Aposs) ksp_Aposs.setType('preonly') pc_Aposs = ksp_Aposs.getPC() pc_Aposs.setType('python') pc_Aposs.setPythonContext(invAposs_ctx(Bs_ksp,projVposs)) ksp_Aposs.setUp() work_1.set(1.) Ms = PETSc.Mat().createPython([works_1.getSizes(), works_1.getSizes()], comm=PETSc.COMM_SELF) Ms.setPythonContext(scaledmats_ctx(Aposs, mus, mus)) Ms.setUp() ksp_Ms = PETSc.KSP().create(comm=PETSc.COMM_SELF) ksp_Ms.setOptionsPrefix("ksp_Ms_") ksp_Ms.setOperators(Ms) ksp_Ms.setType('preonly') pc_Ms = ksp_Ms.getPC() pc_Ms.setType('python') pc_Ms.setPythonContext(scaledmats_ctx(invAposs,invmus,invmus) ) ksp_Ms.setFromOptions() #once a ksp has been passed to SLEPs it cannot be used again so we use a second, identical, ksp for SLEPc as a temporary fix ksp_Ms_forSLEPc = PETSc.KSP().create(comm=PETSc.COMM_SELF) ksp_Ms_forSLEPc.setOptionsPrefix("ksp_Ms_") ksp_Ms_forSLEPc.setOperators(Ms) ksp_Ms_forSLEPc.setType('preonly') pc_Ms_forSLEPc = ksp_Ms_forSLEPc.getPC() pc_Ms_forSLEPc.setType('python') pc_Ms_forSLEPc.setPythonContext(scaledmats_ctx(invAposs,invmus,invmus) ) ksp_Ms_forSLEPc.setFromOptions() # the default local solver is the scaled non assembled local matrix (as in BNN) if self.switchtoASM: Atildes = As if mpi.COMM_WORLD.rank == 0: print('Switch to Additive Schwarz instead of BNN.') ksp_Atildes = PETSc.KSP().create(comm=PETSc.COMM_SELF) ksp_Atildes.setOptionsPrefix("ksp_Atildes_") ksp_Atildes.setOperators(Atildes) ksp_Atildes.setType('preonly') pc_Atildes = ksp_Atildes.getPC() pc_Atildes.setType('cholesky') pc_Atildes.setFactorSolverType('mumps') ksp_Atildes.setFromOptions() #once a ksp has been passed to SLEPs it cannot be used again so we use a second, identical, ksp for SLEPc as a temporary fix ksp_Atildes_forSLEPc = PETSc.KSP().create(comm=PETSc.COMM_SELF) ksp_Atildes_forSLEPc.setOptionsPrefix("ksp_Atildes_") ksp_Atildes_forSLEPc.setOperators(Atildes) ksp_Atildes_forSLEPc.setType('preonly') pc_Atildes_forSLEPc = ksp_Atildes_forSLEPc.getPC() pc_Atildes_forSLEPc.setType('cholesky') pc_Atildes_forSLEPc.setFactorSolverType('mumps') ksp_Atildes_forSLEPc.setFromOptions() if self.switchtoASMpos: if mpi.COMM_WORLD.rank == 0: print('switchtoASMpos has been ignored in favour of switchtoASM.') elif self.switchtoASMpos: Atildes = RsAposRsts if mpi.COMM_WORLD.rank == 0: print('Switch to Apos Additive Schwarz instead of BNN.') ksp_Atildes = PETSc.KSP().create(comm=PETSc.COMM_SELF) ksp_Atildes.setOptionsPrefix("ksp_Atildes_") ksp_Atildes.setOperators(Atildes) ksp_Atildes.setType('preonly') pc_Atildes = ksp_Atildes.getPC() pc_Atildes.setType('python') pc_Atildes.setPythonContext(invRsAposRsts_ctx(As,RsVnegs,RsDnegs,works)) ksp_Atildes.setFromOptions() #once a ksp has been passed to SLEPs it cannot be used again so we use a second, identical, ksp for SLEPc as a temporary fix ksp_Atildes_forSLEPc = PETSc.KSP().create(comm=PETSc.COMM_SELF) ksp_Atildes_forSLEPc.setOptionsPrefix("ksp_Atildes_") ksp_Atildes_forSLEPc.setOperators(Atildes) ksp_Atildes_forSLEPc.setType('preonly') pc_Atildes_forSLEPc = ksp_Atildes_forSLEPc.getPC() pc_Atildes_forSLEPc.setType('python') pc_Atildes_forSLEPc.setPythonContext(invRsAposRsts_ctx(As,RsVnegs,RsDnegs,works)) ksp_Atildes_forSLEPc.setFromOptions() else: #(default) Atildes = Ms ksp_Atildes = PETSc.KSP().create(comm=PETSc.COMM_SELF) ksp_Atildes.setOptionsPrefix("ksp_Atildes_") ksp_Atildes.setOperators(Atildes) ksp_Atildes.setType('preonly') pc_Atildes = ksp_Atildes.getPC() pc_Atildes.setType('python') pc_Atildes.setPythonContext(scaledmats_ctx(invAposs,invmus,invmus) ) ksp_Atildes.setFromOptions() #once a ksp has been passed to SLEPs it cannot be used again so we use a second, identical, ksp for SLEPc as a temporary fix ksp_Atildes_forSLEPc = PETSc.KSP().create(comm=PETSc.COMM_SELF) ksp_Atildes_forSLEPc.setOptionsPrefix("ksp_Atildes_") ksp_Atildes_forSLEPc.setOperators(Atildes) ksp_Atildes_forSLEPc.setType('preonly') pc_Atildes_forSLEPc = ksp_Atildes_forSLEPc.getPC() pc_Atildes_forSLEPc.setType('python') pc_Atildes_forSLEPc.setPythonContext(scaledmats_ctx(invAposs,invmus,invmus) ) ksp_Atildes_forSLEPc.setFromOptions() labs=[] for i, tmp in enumerate(Dnegs): labs.append(f'(\Lambda_-^s)_{i} = {-1.*tmp}') minV0 = minimal_V0(ksp_Atildes,invmusVnegs,labs) #won't compute anything more vecause the solver for Atildes is not mumps minV0s = minV0.V0s labs = minV0.labs if self.viewminV0 == True: minV0.view() self.A = A self.Apos = Apos self.Aneg = Aneg self.Ms = Ms self.As = As self.RsAposRsts = RsAposRsts self.ksp_Atildes = ksp_Atildes self.ksp_Ms = ksp_Ms self.ksp_Atildes_forSLEPc = ksp_Atildes_forSLEPc self.ksp_Ms_forSLEPc = ksp_Ms_forSLEPc self.work = work_1 self.work_2 = work_2 self.works_1 = works_1 self.works_2 = works_2 self.scatter_l2g = scatter_l2g self.mult_max = mult_max self.ksp_Atildes = ksp_Atildes self.minV0 = minV0 self.labs = labs self.Dnegs = Dnegs self.nnegs = nnegs self.works_1.set(1.) self.RsAposRsts.mult(self.works_1,self.works_2) if self.GenEO == True: print(f'{labs=}') self.GenEOV0 = GenEO_V0(self.ksp_Atildes_forSLEPc,self.Ms,self.RsAposRsts,self.mult_max,minV0s,labs,self.ksp_Ms_forSLEPc) self.V0s = self.GenEOV0.V0s if self.viewGenEOV0 == True: self.GenEOV0.view() print(f'{self.GenEOV0.labs=}') else: self.V0s = minV0s self.proj2 = coarse_operators(self.V0s,self.Apos,self.scatter_l2g,self.works_1,self.work) if self.viewV0 == True: self.proj2.view() # work.set(1.) # test = work.copy() # test = self.proj2.coarse_init(work) # testb = work.copy() # self.proj2.project(testb) # testc = work.copy() # self.proj2.project_transpose(testc) # testd = work.copy() # self.apply([], work,testd) self.H2 = PETSc.Mat().createPython([work_1.getSizes(), work_1.getSizes()], comm=PETSc.COMM_WORLD) self.H2.setPythonContext(H2_ctx(self.H2projCS, self.H2addCS, self.proj2, self.scatter_l2g, self.ksp_Atildes, self.works_1, self.works_2 )) self.H2.setUp() self.ksp_Apos = PETSc.KSP().create(comm=PETSc.COMM_WORLD) self.ksp_Apos.setOptionsPrefix("ksp_Apos_") self.ksp_Apos.setOperators(Apos) self.ksp_Apos.setType("cg") if self.compute_ritz_apos: self.ksp_Apos.setComputeEigenvalues(True) self.pc_Apos = self.ksp_Apos.getPC() self.pc_Apos.setType('python') self.pc_Apos.setPythonContext(H2_ctx(self.H2projCS, self.H2addCS, self.proj2, self.scatter_l2g, self.ksp_Atildes, self.works_1, self.works_2 )) self.ksp_Apos.setFromOptions() self.pc_Apos.setFromOptions() #At this point the preconditioner for Apos is ready if self.verbose: if mpi.COMM_WORLD.rank == 0: print(f'#V0(H2) = rank(Ker(Pi2)) = {len(self.proj2.V0)}') Vneg = [] for i in range(mpi.COMM_WORLD.size): nnegi = len(Vnegs) if i == mpi.COMM_WORLD.rank else None nnegi = mpi.COMM_WORLD.bcast(nnegi, root=i) for j in range(nnegi): if i == mpi.COMM_WORLD.rank: works_1 = Vnegs[j].copy() else: works_1.set(0.) self.work.set(0) self.scatter_l2g(works_1, self.work, PETSc.InsertMode.ADD_VALUES) Vneg.append(self.work.copy()) AposinvV0 = [] self.ritz_eigs_apos = None for vec in Vneg: self.ksp_Apos.solve(vec,self.work_2) if self.compute_ritz_apos and self.ritz_eigs_apos is None: self.ritz_eigs_apos = self.ksp_Apos.computeEigenvalues() self.ksp_Apos.setComputeEigenvalues(False) AposinvV0.append(self.work_2.copy()) self.AposinvV0 = AposinvV0 self.proj3 = coarse_operators(self.AposinvV0,self.A,self.scatter_l2g,self.works_1,self.work,V0_is_global=True) self.proj = self.proj3 #this name is consistent with the proj in PCBNN ############################### # ###Alternative to assembling the second coarse operators # # ### # self.Id = PETSc.Mat().createPython([work.getSizes(), work.getSizes()], comm=PETSc.COMM_WORLD) # self.Id.setPythonContext(Id_ctx()) # self.Id.setUp() # # #self.Id = PETSc.Mat().create(comm=PETSc.COMM_SELF) # #self.Id.setType("constantdiagonal") #I don't know how to set the value to 1 # # #self.N = PETSc.Mat().createPython([work.getSizes(), work.getSizes()], comm=PETSc.COMM_WORLD) # #self.N.setPythonContext(N_ctx(self.Aneg,self.A,self.ksp_Apos,self.work,self.work_2)) # #self.N.setUp() # # #self.ksp_N = PETSc.KSP().create(comm=PETSc.COMM_WORLD) # #self.ksp_N.setOptionsPrefix("ksp_N_") # #self.ksp_N.setOperators(self.N) # #self.ksp_N.setType("gmres") # #self.ksp_N.setGMRESRestart(151) ## # if self.compute_ritz_N: # #self.ksp_N.setComputeEigenvalues(True) # ##self.pc_N = self.ksp_N.getPC() # ##self.pc_N.setType('python') # ##self.pc_N.setPythonContext( # #self.ksp_N.setFromOptions() # self.proj4 = coarse_operators(Vneg,self.Id,self.scatter_l2g,self.works_1,self.work,V0_is_global=True) # # self.ProjA = PETSc.Mat().createPython([work.getSizes(), work.getSizes()], comm=PETSc.COMM_WORLD) # self.ProjA.setPythonContext(ProjA_ctx(self.proj4,self.A)) # self.ProjA.setUp() # self.work.set(1.) # #test = self.work.duplicate() # #self.ProjA.mult(self.work,test) # #print('self.ProjA works ok') # # self.ksp_ProjA = PETSc.KSP().create(comm=PETSc.COMM_WORLD) # self.ksp_ProjA.setOptionsPrefix("ksp_ProjA_") # self.ksp_ProjA.setOperators(self.ProjA) # self.ksp_ProjA.setType("gmres") # self.ksp_ProjA.setGMRESRestart(151) # self.ksp_ProjA.setComputeEigenvalues(True) # #self.pc_ProjA = self.ksp_N.getPC() # #self.pc_ProjA.setType('python') # #self.pc_ProjA.setPythonContext( # self.ksp_ProjA.setFromOptions() ############################### ## if self.viewV0 == True: self.proj.view() if self.viewPC == True: self.view() ##Debug DEBUG # works_3 = works.copy() ##projVnegs is a projection # #works.setRandom() # works.set(1.) # projVnegs.mult(works,works_2) # projVnegs.mult(works_2,works_3) # print(f'check that projVnegs is a projection {works_2.norm()} = {works_3.norm()} < {works.norm()}') ##projVposs is a projection ##Pythagoras ok # works.setRandom() # #works.set(1.) # projVnegs.mult(works,works_2) # projVposs.mult(works,works_3) # print(f'{works_2.norm()**2} + {works_3.norm()**2}= {works_2.norm()**2 + works_3.norm()**2} = {(works.norm())**2}') # print(f'0 = {(works - works_2 - works_3).norm()} if the two projections sum to identity') ##Aposs = projVposs Bs projVposs = Bs projVposs (it is implemented as Bs + Anegs) # works_4 = works.copy() # works.setRandom() # #works.set(1.) # projVposs.mult(works,works_2) # Bs.mult(works_2,works_3) # projVposs.mult(works_3,works_2) # Aposs.mult(works,works_4) # print(f'check Aposs = projVposs Bs projVposs = Bs projVposs: {works_2.norm()} = {works_3.norm()} = {works_4.norm()}') # print(f'norms of diffs (should be zero): {(works_2 - works_3).norm()}, {(works_2 - works_4).norm()}, {(works_3 - works_4).norm()}') ###check that Aposs > 0 and Anegs >0 but Bs is indefinite + "Pythagoras" # works_4 = works.copy() # works.set(1.) #(with vector full of ones I get a negative Bs semi-norm) # Bs.mult(works,works_4) # Aposs.mult(works,works_2) # Anegs.mult(works,works_3) # print(f'|.|_Bs {works_4.dot(works)} (can be neg or pos); |.|_Aposs {works_2.dot(works)} > 0; |.|_Anegs {works_3.dot(works)} >0') # print(f' |.|_Bs^2 = |.|_Aposs^2 - |.|_Anegs ^2 = {works_2.dot(works)} - {works_3.dot(works)} = {works_2.dot(works) - works_3.dot(works)} = {works_4.dot(works)} ')## ###check that ksp_Aposs.solve(Aposs * x) = projVposs x # works_4 = works.copy() # works.setRandom() # #works.set(1.) # projVposs.mult(works,works_2) # Aposs(works,works_3) # ksp_Aposs.solve(works_3,works_4) # works_5 = works_2 - works_4 # print(f'norm x = {works.norm()}; norm projVposs x = {works_2.norm()} = norm Aposs\Aposs*x = {works_4.norm()}; normdiff = {works_5.norm()}') ####check that mus*invmus = vec of ones # works.set(1.0) # works_2 = invmus*mus # works_3 = works - works_2 # print(f'0 = norm(vec of ones - mus*invmus) = {works_3.norm()}, mus in [{mus.min()}, {mus.max()}], invmus in [{invmus.min()}, {invmus.max()}]') ###check that Ms*ksp_Ms.solve(Ms*x) = Ms*x # works_4 = works.copy() # works.setRandom() # Atildes.mult(works,works_3) # self.ksp_Atildes.solve(works_3,works_4) # Atildes.mult(works_4,works_2) # works_5 = works_2 - works_3 # print(f'norm x = {works.norm()}; Atilde*x = {works_3.norm()} = norm Atilde*(Atildes\Atildes)*x = {works_2.norm()}; normdiff = {works_5.norm()}') ###check Apos by implementing it a different way in Apos_debug # Apos_debug = PETSc.Mat().createPython([work.getSizes(), work.getSizes()], comm=PETSc.COMM_WORLD) # Apos_debug.setPythonContext(Apos_debug_ctx(projVposs, Aposs, scatter_l2g, works, work)) # Apos_debug.setUp() # work.setRandom() # test = work.duplicate() # test2 = work.duplicate() # Apos.mult(work,test) # Apos_debug.mult(work,test2) # testdiff = test-test2 # print(f'norm of |.|_Apos = {np.sqrt(test.dot(work))} = |.|_Apos_debug = {np.sqrt(test2.dot(work))} ; norm of diff = {testdiff.norm()}') ### ###check that the projection in proj2 is a self.proj2.A orth projection #work.setRandom() # work.set(1.) # test = work.copy() # self.proj2.project(test) # test2 = test.copy() # self.proj2.project(test2) # testdiff = test-test2 # print(f'norm(Pi x - Pi Pix) = {testdiff.norm()} = 0') # self.proj2.A.mult(test,test2) # test3 = work.duplicate() # self.proj2.A.mult(work,test3) # print(f'|Pi x|_A^2 - |x|_A^2 = {test.dot(test2)} - {work.dot(test3)} = {test.dot(test2) - work.dot(test3)} < 0 ') # #test2 = A Pi x ( = Pit A Pi x) # test3 = test2.copy() # self.proj2.project_transpose(test3) # test = test3.copy() # self.proj2.project_transpose(test) # testdiff = test3 - test2 # print(f'norm(A Pi x - Pit A Pix) = {testdiff.norm()} = 0 = {(test - test3).norm()} = norm(Pit Pit A Pi x - Pit A Pix); compare to norm(A Pi x) = {test2.norm()} ') # #work.setRandom() # work.set(1.) # test2 = work.copy() # self.proj2.project_transpose(test2) # test2 = -1*test2 # test2 += work # # test = work.copy() # test = self.proj2.coarse_init(work) # test3 = work.duplicate() # self.proj2.A.mult(test,test3) ###check that the projection in proj3 is a self.proj3.A orth projection whose image includes Ker(Aneg) # #work.setRandom() # work.set(1.) # test = work.copy() # self.proj3.project(test) # test2 = test.copy() # self.proj3.project(test2) # testdiff = test-test2 # print(f'norm(Pi x - Pi Pix) = {testdiff.norm()} = 0') # self.proj3.A.mult(test,test2) # test3 = work.duplicate() # self.proj3.A.mult(work,test3) # print(f'|Pi x|_A^2 - |x|_A^2 = {test.dot(test2)} - {work.dot(test3)} = {test.dot(test2) - work.dot(test3)} < 0 ') # #test2 = A Pi x ( = Pit A Pi x) # test3 = test2.copy() # self.proj3.project_transpose(test3) # test = test3.copy() # self.proj3.project_transpose(test) # testdiff = test3 - test2 # print(f'norm(A Pi x - Pit A Pix) = {testdiff.norm()} = 0 = {(test - test3).norm()} = norm(Pit Pit A Pi x - Pit A Pix); compare to norm(A Pi x) = {test2.norm()} ') # #work.setRandom() # work.set(1.) # test2 = work.copy() # self.proj3.project_transpose(test2) # test2 = -1*test2 # test2 += work # # test = work.copy() # test = self.proj3.coarse_init(work) # test3 = work.duplicate() # self.proj3.A.mult(test,test3) # # print(f'norm(A coarse_init(b)) = {test3.norm()} = {test2.norm()} = norm((I-Pit b)); norm diff = {(test2 - test3).norm()}') # # work.set(1.) # test = work.copy() # test2 = work.copy() # self.proj3.project(test2) # test3 = work.copy() # self.proj3.project(test3) # test = work.copy() # self.Apos.mult(test2,test) # test2 = work.copy() # self.A.mult(test3,test2) # print(f'norm(Apos Pi3 x) = {test.norm()} = {test2.norm()} = norm(A Pi3 x); norm diff = {(test - test2).norm()}') # for vec in self.AposinvV0: # test = vec.copy() # self.proj3.project(test) # print(f'norm(Pi3 AposinvV0[i]) = {test.norm()} compare to norm of the non projected vector norm ={(vec).norm()}') # ### END Debug DEBUG def mult(self, x, y): """ Applies the domain decomposition preconditioner followed by the projection preconditioner to a vector. Parameters ========== x : petsc.Vec The vector to which the preconditioner is to be applied. y : petsc.Vec The vector that stores the result of the preconditioning operation. """ ######################## ######################## xd = x.copy() if self.H3projCS == True: self.proj3.project_transpose(xd) self.H2.mult(xd,y) if self.H3projCS == True: self.proj3.project(y) if self.H3addCS == True: xd = x.copy() ytild = self.proj3.coarse_init(xd) # I could save a coarse solve by combining this line with project_transpose if ytild.dot(xd) < 0: print(f'x.dot(coarse_init(x)) = {ytild.dot(xd)} < 0 ') y += ytild def MP_mult(self, x, y): """ Applies the domain decomposition multipreconditioner followed by the projection preconditioner to a vector. Parameters ========== x : petsc.Vec The vector to which the preconditioner is to be applied. y : FIX The list of ndom vectors that stores the result of the multipreconditioning operation (one vector per subdomain). """ print('not implemented') def apply(self, pc, x, y): """ Applies the domain decomposition preconditioner followed by the projection preconditioner to a vector. This is just a call to PCNew.mult with the function name and arguments that allow PCNew to be passed as a preconditioner to PETSc.ksp. Parameters ========== pc: This argument is not called within the function but it belongs to the standard way of calling a preconditioner. x : petsc.Vec The vector to which the preconditioner is to be applied. y : petsc.Vec The vector that stores the result of the preconditioning operation. """ self.mult(x,y) def view(self): self.gathered_ns = mpi.COMM_WORLD.gather(self.ns, root=0) self.gathered_nints = mpi.COMM_WORLD.gather(self.nints, root=0) self.gathered_Gammas = mpi.COMM_WORLD.gather(self.nGammas, root=0) self.minV0.gathered_dim = mpi.COMM_WORLD.gather(self.minV0.nrb, root=0) self.gathered_labs = mpi.COMM_WORLD.gather(self.labs, root=0) self.gathered_nneg = mpi.COMM_WORLD.gather(self.nnegs, root=0) self.gathered_Dneg = mpi.COMM_WORLD.gather(self.Dnegs, root=0) if self.GenEO == True: self.GenEOV0.gathered_nsharp = mpi.COMM_WORLD.gather(self.GenEOV0.n_GenEO_eigmax, root=0) self.GenEOV0.gathered_nflat = mpi.COMM_WORLD.gather(self.GenEOV0.n_GenEO_eigmin, root=0) self.GenEOV0.gathered_dimKerMs = mpi.COMM_WORLD.gather(self.GenEOV0.dimKerMs, root=0) self.GenEOV0.gathered_Lambdasharp = mpi.COMM_WORLD.gather(self.GenEOV0.Lambda_GenEO_eigmax, root=0) self.GenEOV0.gathered_Lambdaflat = mpi.COMM_WORLD.gather(self.GenEOV0.Lambda_GenEO_eigmin, root=0) if mpi.COMM_WORLD.rank == 0: print('#############################') print(f'view of PCNew') print(f'{self.switchtoASM=}') print(f'{self.verbose= }') print(f'{self.GenEO= }') print(f'{self.H3addCS= }') print(f'{self.H3projCS= }') print(f'{self.H2projCS= }') print(f'{self.viewPC= }') print(f'{self.viewV0= }') print(f'{self.viewGenEOV0= }') print(f'{self.viewnegV0= }') print(f'{self.viewminV0= }') print(f'{self.compute_ritz_apos=}') print(f'{self.mult_max=}') print(f'### info about the subdomains ###') self.nint = np.sum(self.gathered_nints) self.nGamma = self.nglob - self.nint print(f'{self.gathered_ns =}') print(f'{self.gathered_nints =}') print(f'{self.gathered_Gammas=}') print(f'{self.nGamma=}') print(f'{self.nint=}') print(f'{self.nglob=}') print(f'{self.gathered_labs=}') print(f'### info about minV0.V0s = (Ker(Atildes)) ###') print(f'{self.minV0.mumpsCntl3=}') print(f'###info about Vnegs = rank(Anegs) = coarse components for proj3') print(f'{self.gathered_nneg=}') print(f'{np.sum(self.gathered_nneg)=}') if (self.ksp_Atildes.pc.getFactorSolverType() == 'mumps'): print(f'dim(Ker(Atildes)) = {self.minV0.gathered_dim}') else: print(f'Ker(Atildes) not computed because pc is not mumps') if self.GenEO == True: print(f'### info about GenEOV0.V0s ###') print(f'{self.GenEOV0.tau_eigmax=}') print(f'{self.GenEOV0.tau_eigmin=}') print(f'{self.GenEOV0.eigmax=}') print(f'{self.GenEOV0.eigmin=}') print(f'{self.GenEOV0.nev=}') print(f'{self.GenEOV0.maxev=}') print(f'{self.GenEOV0.mumpsCntl3=}') print(f'{self.GenEOV0.verbose=}') print(f'{self.GenEOV0.gathered_nsharp=}') print(f'{self.GenEOV0.gathered_nflat=}') #print(f'{self.GenEOV0.gathered_dimKerMs=}') #print(f'{np.array(self.GenEOV0.gathered_Lambdasharp)=}') #print(f'{np.array(self.GenEOV0.gathered_Lambdaflat)=}') print(f'### info about the preconditioner for Apos ###') print(f'{self.proj2.V0_is_global=}') if(self.proj2.V0_is_global == False): print(f'{self.proj2.gathered_dimV0s=}') if self.GenEO == True: print(f'global dim V0 for Apos = {self.proj2.dim} = ({np.sum(self.gathered_nneg)} from Vneg ) + ({np.sum(self.minV0.gathered_dim)} from Ker(Atildes)) + ({np.sum(self.GenEOV0.gathered_nsharp)} from GenEO_eigmax) + ({np.sum(self.GenEOV0.gathered_nflat) } from GenEO_eigmin)') else: print(f'global dim V0 for Apos = {np.sum(self.proj2.gathered_dimV0s)} = ({np.sum(self.minV0.gathered_dim)} from Ker(Atildes))') if self.compute_ritz_apos and self.ritz_eigs_apos is not None: print(f'Estimated kappa(H2 Apos) = {self.ritz_eigs_apos.max()/self.ritz_eigs_apos.min() }; with lambdamin = {self.ritz_eigs_apos.min()} and lambdamax = {self.ritz_eigs_apos.max()}') print('#############################') self.savetofile() def savetofile(self): if mpi.COMM_WORLD.rank == 0: if not os.path.exists(self.test_case): os.mkdir(self.test_case) np.savez(f'{self.test_case}/init', switchtoASM = self.switchtoASM, verbose = self.verbose, GenEO = self.GenEO, H3addCS = self.H3addCS, H3projCS = self.H3projCS, H2projCS = self.H2projCS, viewPC = self.viewPC, viewV0 = self.viewV0, viewGenEOV0 = self.viewGenEOV0, viewnegV0 = self.viewnegV0, viewminV0 = self.viewminV0, compute_ritz_apos = self.compute_ritz_apos, mult_max = self.mult_max, gathered_ns = self.gathered_ns, gathered_nints = self.gathered_nints, gathered_Gammas = self.gathered_Gammas, nGamma = self.nGamma, nint = self.nint, nglob = self.nglob, minV0_mumpsCntl3 = self.minV0.mumpsCntl3, gathered_labs= np.asarray(self.gathered_labs,dtype='object'), gathered_nneg = self.gathered_nneg, minV0_gathered_dim = self.minV0.gathered_dim, ritz_eigs_Apos = self.ritz_eigs_apos , sum_nneg = np.sum(self.gathered_nneg), proj2_V0_is_global = self.proj2.V0_is_global, proj2_gathered_dimV0s = np.asarray(self.proj2.gathered_dimV0s), proj2_dimV0 = np.sum(self.proj2.gathered_dimV0s), proj2_sum_dimminV0 = np.sum(self.minV0.gathered_dim) , ) if self.GenEO == True: np.savez(f'{self.test_case}/GenEO', GenEOV0_tau_eigmax = self.GenEOV0.tau_eigmax, GenEOV0_tau_eigmin = self.GenEOV0.tau_eigmin, GenEOV0_eigmax = self.GenEOV0.eigmax, GenEOV0_eigmin = self.GenEOV0.eigmin, GenEOV0_nev = self.GenEOV0.nev, GenEOV0_maxev = self.GenEOV0.maxev, GenEOV0_mumpsCntl3 = self.GenEOV0.mumpsCntl3, GenEOV0_verbose = self.GenEOV0.verbose, GenEOV0_gathered_nsharp = np.asarray(self.GenEOV0.gathered_nsharp), GenEOV0_gathered_nflat = np.asarray(self.GenEOV0.gathered_nflat), GenEOV0_sum_nsharp = np.sum(self.GenEOV0.gathered_nsharp), GenEOV0_sum_nflat = np.sum(self.GenEOV0.gathered_nflat), ) class Aneg_ctx(object): def __init__(self, Vnegs, DVnegs, scatter_l2g, works, works_2): self.scatter_l2g = scatter_l2g self.works = works self.works_2 = works_2 self.Vnegs = Vnegs self.DVnegs = DVnegs self.gamma = PETSc.Vec().create(comm=PETSc.COMM_SELF) self.gamma.setType(PETSc.Vec.Type.SEQ) self.gamma.setSizes(len(self.Vnegs)) def mult(self, mat, x, y): y.set(0) self.scatter_l2g(x, self.works, PETSc.InsertMode.INSERT_VALUES, PETSc.ScatterMode.SCATTER_REVERSE) self.works_2.set(0) for i, vec in enumerate(self.Vnegs): self.works_2.axpy(self.works.dot(self.DVnegs[i]), vec) self.scatter_l2g(self.works_2, y, PETSc.InsertMode.ADD_VALUES) class Apos_debug_ctx(object): def __init__(self, projVposs, Aposs, scatter_l2g, works, work): self.scatter_l2g = scatter_l2g self.work = works self.works = works self.projVposs = projVposs self.Aposs = Aposs def mult(self, mat, x, y): y.set(0) works_2 = self.works.duplicate() self.scatter_l2g(x, self.works, PETSc.InsertMode.INSERT_VALUES, PETSc.ScatterMode.SCATTER_REVERSE) self.Aposs.mult(self.works,works_2) self.scatter_l2g(works_2, y, PETSc.InsertMode.ADD_VALUES) class Apos_ctx(object): def __init__(self,A, Aneg): self.A = A self.Aneg = Aneg def mult(self, mat, x, y): xtemp = x.duplicate() self.Aneg.mult(x,xtemp) self.A.mult(x,y) y += xtemp class Anegs_ctx(object): def __init__(self, Vnegs, DVnegs): self.Vnegs = Vnegs self.DVnegs = DVnegs def mult(self, mat, x, y): y.set(0) for i,vec in enumerate(self.Vnegs): y.axpy(x.dot(self.DVnegs[i]), vec) class RsAposRsts_ctx(object): def __init__(self,As,RsVnegs,RsDVnegs): self.As = As self.RsVnegs = RsVnegs self.RsDVnegs = RsDVnegs def mult(self, mat, x, y): self.As.mult(x,y) for i,vec in enumerate(self.RsVnegs): y.axpy(x.dot(self.RsDVnegs[i]), vec) class invRsAposRsts_ctx(object): def __init__(self,As,RsVnegs,RsDnegs,works): self.As = As self.works = works self.RsVnegs = RsVnegs self.RsDnegs = RsDnegs self.ksp_As = PETSc.KSP().create(comm=PETSc.COMM_SELF) self.ksp_As.setOptionsPrefix("ksp_As_") self.ksp_As.setOperators(self.As) self.ksp_As.setType('preonly') self.pc_As = self.ksp_As.getPC() self.pc_As.setType('cholesky') self.pc_As.setFactorSolverType('mumps') self.ksp_As.setFromOptions() self.AsinvRsVnegs = [] for i,vec in enumerate(self.RsVnegs): self.ksp_As.solve(vec,self.works) self.AsinvRsVnegs.append(self.works.copy()) self.Matwood = PETSc.Mat().create(comm=PETSc.COMM_SELF) self.Matwood.setType(PETSc.Mat.Type.SEQDENSE) self.Matwood.setSizes([len(self.AsinvRsVnegs),len(self.AsinvRsVnegs)]) self.Matwood.setOption(PETSc.Mat.Option.SYMMETRIC, True) self.Matwood.setPreallocationDense(None) for i, vec in enumerate(self.AsinvRsVnegs): for j in range(i): tmp = self.RsVnegs[j].dot(vec) self.Matwood[i, j] = tmp self.Matwood[j, i] = tmp tmp = self.RsVnegs[i].dot(vec) self.Matwood[i, i] = tmp + 1/self.RsDnegs[i] self.Matwood.assemble() self.ksp_Matwood = PETSc.KSP().create(comm=PETSc.COMM_SELF) self.ksp_Matwood.setOperators(self.Matwood) self.ksp_Matwood.setType('preonly') self.pc = self.ksp_Matwood.getPC() self.pc.setType('cholesky') self.gamma, _ = self.Matwood.getVecs() self.alpha = self.gamma.duplicate() def mult(self, mat, x, y): self.ksp_As.solve(x,y) for i, vec in enumerate(self.AsinvRsVnegs): self.gamma[i] = vec.dot(x) self.ksp_Matwood.solve(self.gamma, self.alpha) for i, vec in enumerate(self.AsinvRsVnegs): y.axpy(-self.alpha[i], vec) def apply(self,pc, x, y): self.mult(pc,x,y) class Aposs_ctx(object): def __init__(self,Bs, Anegs): self.Bs = Bs self.Anegs = Anegs def mult(self, mat, x, y): xtemp = x.duplicate() self.Anegs.mult(x,xtemp) self.Bs.mult(x,y) y += xtemp class scaledmats_ctx(object): def __init__(self, mats, musl, musr): self.mats = mats self.musl = musl self.musr = musr def mult(self, mat, x, y): xtemp = x.copy()*self.musr self.mats.mult(xtemp,y) y *= self.musl def apply(self, mat, x, y): self.mult(mat, x, y) class invAposs_ctx(object): def __init__(self,Bs_ksp,projVposs): self.Bs_ksp = Bs_ksp self.projVposs = projVposs def apply(self, mat, x, y): xtemp1 = y.duplicate() xtemp2 = y.duplicate() self.projVposs.mult(x,xtemp1) self.Bs_ksp.solve(xtemp1,xtemp2) self.projVposs.mult(xtemp2,y) def mult(self, mat, x, y): #xtemp1 = y.duplicate() #xtemp2 = y.duplicate() #self.projVnegs.mult(x,xtemp1) #self.Bs_ksp.solve(xtemp1,xtemp2) #self.projVnegs.mult(xtemp2,y) self.apply(mat, x, y) class projVnegs_ctx(object): def __init__(self, Vnegs): self.Vnegs = Vnegs def mult(self, mat, x, y): y.set(0) for i,vec in enumerate(self.Vnegs): y.axpy(x.dot(vec) , vec) class projVposs_ctx(object): def __init__(self, projVnegs): self.projVnegs = projVnegs def mult(self, mat, x, y): self.projVnegs(-x,y) y.axpy(1.,x) class H2_ctx(object): def __init__(self, projCS, addCS, proj2, scatter_l2g, ksp_Atildes, works_1, works_2 ): self.projCS = projCS self.addCS = addCS self.proj2 = proj2 self.scatter_l2g = scatter_l2g self.ksp_Atildes = ksp_Atildes self.works_1 = works_1 self.works_2 = works_2 def mult(self,mat,x,y): self.apply([],x,y) def apply(self,pc, x, y): xd = x.copy() if self.projCS == True: self.proj2.project_transpose(xd) self.scatter_l2g(xd, self.works_1, PETSc.InsertMode.INSERT_VALUES, PETSc.ScatterMode.SCATTER_REVERSE) self.ksp_Atildes.solve(self.works_1, self.works_2) y.set(0.) self.scatter_l2g(self.works_2, y, PETSc.InsertMode.ADD_VALUES) if self.projCS == True: self.proj2.project(y) if self.addCS == True: xd = x.copy() ytild = self.proj2.coarse_init(xd) # I could save a coarse solve by combining this line with project_transpose #print(f'in H2 x.dot(coarse_init(x)) = {ytild.dot(xd)} > 0 ') if ytild.dot(xd) < 0: print('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') y += ytild class N_ctx(object): #(I - Aneg *invApos) def __init__(self,Aneg,A,ksp_Apos,work,work_2): self.Aneg = Aneg self.A = A self.ksp_Apos = ksp_Apos self.work = work self.work_2 = work_2 def mult(self, mat, x, y): if mpi.COMM_WORLD.rank == 0: print('in N_ctx') self.ksp_Apos.solve(x,self.work) #self.A.mult(self.work,y) self.Aneg.mult(self.work,self.work_2) #self.Aneg.mult(self.work,y) y.set(0.) y.axpy(1.,x) y.axpy(-1.,self.work_2) class Id_ctx(object): # I # def __init__(self): def mult(self, mat, x, y): y.axpy(1.,x) class ProjA_ctx(object): #(Vneg (Vnegt *Vneg)\Vnegt *A ) def __init__(self, proj4, A): self.proj4 = proj4 self.A = A def mult(self, mat, x, y): xd = x.copy() self.proj4.project(xd) self.A.mult(xd,y) ###UNFINISHED #class Psi_ctx(object) #Dneg^{-1} - Vnegt Aneg^{-1} Vneg # def __init__(self,Vneg,ksp_Apos,work,work_2): # self.Vneg = Vneg # self.work = work # self.work = work_2 # #self.gamma = PETSc.Vec().create(comm=PETSc.COMM_SELF) # #self.gamma.setType(PETSc.Vec.Type.SEQ) # #self.gamma.setSizes(len(self.Vneg)) # self.ksp_Apos = ksp_Apos # def mult(self, mat, x, y): # part with Dneg inv is not at all implemented yet # self.work.set(0.) # for i, vec in enumerate(self.Vneg): # self.work.axpy(x[i], vec) # # self.ksp_Apos.solve(self.work,self.work_2) # # #y = self.gamma.duplicate() # for i, vec in enumerate(self.V0): # y[i] = vec.dot(x) #