It appears there are currently no Python bindings for > https://petsc.org/main/manualpages/Mat/MatSetPreallocationCOOLocal/ > https://petsc.org/main/manualpages/Mat/MatSetValuesCOO/
they should be fairly easy to add in a merge request. Barry > On Aug 17, 2023, at 5:57 PM, Matthew Knepley <[email protected]> wrote: > > On Fri, Aug 18, 2023 at 12:49 AM Erik Kneller via petsc-users > <[email protected] <mailto:[email protected]>> wrote: >> Hi All, >> >> I need to fill non-zero values of a Petsc matrix via petsc4py for the domain >> defined by A.getOwnershipRange() using three Numpy arrays: (1) array >> containing row indices of non-zero value, (2) array containing column >> indices of non-zero values and (3) array containing the non-zero matrix >> values. How can one perform this type of filling operation in petsc4py? The >> method A.setValues does not appear to allow this since it only works on an >> individual matrix element or a block of matrix elements. >> >> I am using Numpy arrays since they can be computed in loops optimized using >> Numba on each processor. I also cannot pass the Petsc matrix to a Numba >> compiled function since type information cannot be inferred. I absolutely >> need to avoid looping in standard Python to define Petsc matrix elements due >> to performance issues. I also need to use a standard petscy4py method and >> avoid writing new C or Fortran wrappers to minimize language complexity. >> >> Example Algorithm Building on Lisandro Dalcin's 2D Poisson Example: >> ---------------------------------------------------------------------------------------------- >> comm = PETSc.COMM_WORLD >> rank = comm.getRank() >> >> dx = 1.0/(xnodes + 1) # xnodes is the number of nodes in the x and >> y-directions of the grid >> nnz_max = 5 # max number of non-zero values per row >> >> A = PETSc.Mat() >> A.create(comm=PETSc.COMM_WORLD) >> A.setSizes((xnodes*ynodes, xnodes*ynodes)) >> A.setType(PETSc.Mat.Type.AIJ) >> A.setPreallocationNNZ(nnz_max) >> >> rstart, rend = A.getOwnershipRange() >> >> # Here Anz, Arow and Acol are vectors with size equal to the number of >> non-zero values >> Anz, Arow, Acol = build_nonzero_numpy_arrays_using_numba(rstart, rend, >> nnz_max, dx, xnodes, ynodes) >> >> A.setValues(Arow, Acol, Anz) # <--- This does not work. > > I see at least two options > > https://petsc.org/main/manualpages/Mat/MatCreateSeqAIJWithArrays/ > > or > > https://petsc.org/main/manualpages/Mat/MatSetPreallocationCOOLocal/ > https://petsc.org/main/manualpages/Mat/MatSetValuesCOO/ > > Thanks, > > Matt > >> A.assemblyBegin() >> A.assemblyEnd() >> ksp = PETSc.KSP() >> ksp.create(comm=A.getComm()) >> ksp.setType(PETSc.KSP.Type.CG <http://petsc.ksp.type.cg/>) >> ksp.getPC().setType(PETSc.PC.Type.GAMG) >> ksp.setOperators(A) >> ksp.setFromOptions() >> x, b = A.createVecs() >> b.set(1.0) >> >> ksp.solve(b, x) >> >> Regards, >> Erik Kneller, Ph.D. >> > > > -- > What most experimenters take for granted before they begin their experiments > is infinitely more interesting than any results to which their experiments > lead. > -- Norbert Wiener > > https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
