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