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
> <[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.
>
> 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.
>
