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

[Zhang, Hong]

You can use this to build a PETSc matrix with the index arrays ai,aj and the value array aa:
    PETSc.Mat().createAIJ(size=(nrows,ncols), csr=(ai,aj,aa))

Hong (Mr.)

> On Aug 17, 2023, at 7:37 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.
> 
> 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.
>