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

[Barry Smith]

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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20230817/89056013/attachment.html>