On Fri, Aug 18, 2023 at 12:49 AM Erik Kneller via petsc-users < petsc-users@mcs.anl.gov> 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) > 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/>