<div dir="ltr"><div dir="ltr">On Fri, Aug 18, 2023 at 12:49 AM Erik Kneller via petsc-users <<a href="mailto:petsc-users@mcs.anl.gov">petsc-users@mcs.anl.gov</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div><div style="font-family:"Helvetica Neue",Helvetica,Arial,sans-serif;font-size:13px"><div dir="ltr">Hi All,</div><div dir="ltr"><br></div><div dir="ltr">I need to fill non-zero values of a Petsc matrix via petsc4py for the domain defined by <span>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 <span>petsc4py</span>? 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.<br></span></div><div dir="ltr"><span><br></span></div><div dir="ltr"><span>I am using Numpy arrays since they can be computed in loops optimized using Numba on each processor. <span><span>I also cannot pass the Petsc matrix to a Numba compiled function since type information cannot be inferred.</span></span> 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.<br></span></div><div dir="ltr"><span><br></span></div><div dir="ltr"><span>Example Algorithm Building on <span>Lisandro Dalcin's 2D Poisson Example:</span></span></div><div dir="ltr"><span><span>----------------------------------------------------------------------------------------------<br></span></span></div><div dir="ltr"><span><span>comm = PETSc.COMM_WORLD</span><br></span></div><div dir="ltr"><span><span>rank = comm.getRank()</span><br></span></div><div dir="ltr"><span><span><br></span></span></div><div dir="ltr"><span><span>dx = 1.0/(xnodes + 1)</span> # xnodes is the number of nodes in the x and y-directions of the grid<br></span></div><div dir="ltr"><span><span><span>nnz_max</span></span> = 5 # max number of non-zero values per row<br></span></div><div dir="ltr"><span><span><br></span></span></div><div dir="ltr"><span><span>A = PETSc.Mat()</span><br></span></div><div dir="ltr"><div>A.create(comm=PETSc.COMM_WORLD)<br>A.setSizes((xnodes*ynodes, xnodes*ynodes))<br><div>A.setType(PETSc.Mat.Type.AIJ)</div><div dir="ltr"><span>A.setPreallocationNNZ(<span><span><span><span>nnz_max</span></span></span></span>)</span></div><div dir="ltr"><span><span><br></span></span></div><div dir="ltr"><span><span>rstart, rend = A.getOwnershipRange()</span></span></div><div dir="ltr"><span><div><div><br></div><div dir="ltr"># Here Anz, Arow and Acol are vectors with size equal to the number of non-zero values <br></div><div>Anz, Arow, Acol = build_nonzero_numpy_arrays_using_numba(rstart, rend, nnz_max, dx, xnodes, ynodes)</div><div><br></div><div>A.setValues(Arow, Acol, Anz) # <--- This does not work.</div></div></span></div></div></div></div></div></blockquote><div><br></div><div>I see at least two options<br><br>  <a href="https://petsc.org/main/manualpages/Mat/MatCreateSeqAIJWithArrays/">https://petsc.org/main/manualpages/Mat/MatCreateSeqAIJWithArrays/</a></div><div><br></div><div>or</div><div><br></div>  <a href="https://petsc.org/main/manualpages/Mat/MatSetPreallocationCOOLocal/">https://petsc.org/main/manualpages/Mat/MatSetPreallocationCOOLocal/</a><br><div>  <a href="https://petsc.org/main/manualpages/Mat/MatSetValuesCOO/">https://petsc.org/main/manualpages/Mat/MatSetValuesCOO/</a></div><div><br></div><div>  Thanks,</div><div><br></div><div>     Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div><div style="font-family:"Helvetica Neue",Helvetica,Arial,sans-serif;font-size:13px"><div dir="ltr"><div><div dir="ltr"><div dir="ltr"><div>A.assemblyBegin()<br><div>A.assemblyEnd()</div><div dir="ltr"><div>ksp = PETSc.KSP()<br>ksp.create(comm=A.getComm())</div><div dir="ltr"><div>ksp.setType(<a href="http://PETSc.KSP.Type.CG" target="_blank">PETSc.KSP.Type.CG</a>)<br>ksp.getPC().setType(PETSc.PC.Type.GAMG)</div><div dir="ltr"><div>ksp.setOperators(A)<br>ksp.setFromOptions()</div><div dir="ltr"><div>x, b = A.createVecs()<br><div>b.set(1.0)</div><div><br></div></div><div dir="ltr"><span>ksp.solve(b, x)</span></div></div></div></div><div><br></div><div dir="ltr">Regards,</div><div dir="ltr">Erik Kneller, Ph.D.<br></div></div></div><div><br></div></div><span></span></div></div></div></div></div></blockquote></div><br clear="all"><div><br></div><span class="gmail_signature_prefix">-- </span><br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>