[petsc-users] Filling non-zero values of a Petsc matrix using numpy arrays with non-zero indices and values (petsc4py)

[Matthew Knepley]

On Fri, Aug 18, 2023 at 12:49 AM Erik Kneller via petsc-users <
petsc-users at 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/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20230818/36bc4eaf/attachment.html>