<div dir="ltr"><div dir="ltr">On Sun, Aug 20, 2023 at 12:53 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></div>
<div dir="ltr">Hi All,</div><div dir="ltr"><br></div><div dir="ltr">Thank you for the recommendations. The first option provided by Matthew works in petsc4py but required some additional computations to re-organize information. I developed an approach that was easier to build within my current system using the option provided by Hong along with some other posts I found on the user forum. See below for an example of how I integrated Numba and petcs4py so the filling of matrix elements on each processor is done using the optimized machine code produced by Numba (you can eliminate the cleaning step if the number of non-zero elements for each row is well defined). Scaling and overall performance is now satisfactory. <br></div><div dir="ltr"><br></div><div dir="ltr">I do have another question. In order to make this work I create the Petsc matrix 'A' twice, first to get local ownership and second to define the local elements for a particular processor:<br></div><div dir="ltr"><br></div><div dir="ltr">Algorithm</div><div dir="ltr">------------<br></div><div dir="ltr">(1) Create a Petsc matrix 'A' and set size and type</div><div dir="ltr">(2) Get the local row start and end for matrix 'A'</div><div dir="ltr">(3) Define the local non-zero coefficients for the rows owned by processor using a Numba JIT-compiled loop and store result in a csr matrix defined using Scipy. <br></div><div dir="ltr">(4) Redefine Petsc matix 'A' using the the local csr elements define in step 3.</div><div dir="ltr">(5) Begin and end assembly</div><div dir="ltr">(6) Define RHS and initial solution vectors</div><div dir="ltr">(7) Solve system <br></div><div dir="ltr"><span> (see code example below for details)</span></div><div dir="ltr"><span><br></span></div><div dir="ltr"><span>Steps (1) and (4) appear redundant and potentially sub-optimal from a performance perspective (perhaps due to my limited experience with petscy4py). Is there a better approach in terms of elegance and performance? </span></div></div></div></blockquote><div><br></div><div>In 1), if you just want the default ownership, you can call</div><div><br></div><div> <a href="https://petsc.org/main/manualpages/Sys/PetscSplitOwnership/">https://petsc.org/main/manualpages/Sys/PetscSplitOwnership/</a></div><div><br></div><div>which is what the Mat is calling underneath.</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"><span>Also, if this is the optimal syntax could someone elaborate on what exactly is occurring behind the seen in terms of memory allocation?<br></span></div><div dir="ltr"><span><br></span></div><div dir="ltr"><span>Thank you all again for your guidance and insights.<br></span></div><div dir="ltr"><span><br></span></div><div dir="ltr"><span><span><span><span>Modified Version of Dalcin's 2D Poisson Example</span></span></span> with Numba<br></span></div><div dir="ltr">----------------------------------------------------------------------------------<br></div><div dir="ltr"><div dir="ltr">import sys<br>from typing import Tuple<br>import numpy.typing as npt<br>import numpy as np<br>from numba import njit<br>import petsc4py<br>petsc4py.init(sys.argv)<br>from petsc4py import PETSc<br>import scipy.sparse as sps<br><br><br>def main() -> None:<br> xnodes = 8000<br> nnz_max = 10<br> <br> ynodes = xnodes<br> N = xnodes*ynodes<br> dx = 1.0/(xnodes + 1)<br> <br> A = PETSc.Mat()<br> A.create(comm=PETSc.COMM_WORLD)<br> A.setSizes((xnodes*ynodes, xnodes*ynodes))<br> A.setType(PETSc.Mat.Type.AIJ)<br> <br> rstart, rend = A.getOwnershipRange()<br> <div> Ascipy = build_csr_matrix(<br> rstart, rend, nnz_max, dx, xnodes, ynodes)<br> csr=(<br> Ascipy.indptr[rstart:rend+1] - Ascipy.indptr[rstart],<br> Ascipy.indices[Ascipy.indptr[rstart]:Ascipy.indptr[rend]],<br> Ascipy.data[Ascipy.indptr[rstart]:Ascipy.indptr[rend]]<br> )<br> A = PETSc.Mat().createAIJ(size=(N,N), csr=csr)</div> <br> A.assemblyBegin()<br> A.assemblyEnd()<br> <br> ksp = PETSc.KSP()<br> ksp.create(comm=A.getComm())<br> ksp.setType(<a href="http://PETSc.KSP.Type.CG" target="_blank">PETSc.KSP.Type.CG</a>)<br> ksp.getPC().setType(PETSc.PC.Type.GAMG)<br> <br> ksp.setOperators(A)<br> ksp.setFromOptions()<br> <br> x, b = A.createVecs()<br> b.set(1.0)<br> <br> ksp.solve(b, x)<br> <br> <br>def build_csr_matrix(<br> rstart:int, <br> rend:int, <br> nnz_max:int, <br> dx:float, <br> xnodes:int, <br> ynodes:int<br>):<br> Anz, Arow, Acol = build_nonzero_arrays(<br> rstart, rend, nnz_max, dx, xnodes, ynodes<br> )<br> N = xnodes*ynodes<br> Ls = sps.csr_matrix((Anz, (Arow,Acol)), shape=(N,N), dtype=np.float64)<br> return Ls<br> <br> <br>def build_nonzero_arrays(<br> rstart:int, <br> rend:int, <br> nnz_max:int, <br> dx:float, <br> xnodes:int, <br> ynodes:int<br>) -> Tuple[<br> npt.NDArray[np.float64], <br> npt.NDArray[np.float64], <br> npt.NDArray[np.float64]<br> ]:<br> nrows_local = (rend - rstart) + 1<br> Anz_ini = np.zeros((nnz_max*nrows_local), dtype=np.float64)<br> Arow_ini = np.zeros((nnz_max*nrows_local), dtype=np.int32)<br> Acol_ini = np.zeros((nnz_max*nrows_local), dtype=np.int32)<br> icount_nz = define_nonzero_values(<br> rstart, rend, dx, xnodes, ynodes, Anz_ini, Arow_ini, Acol_ini<br> )<br> (<br> Anz, Arow, Acol<br> ) = clean_nonzero_arrays(icount_nz, Anz_ini, Arow_ini, Acol_ini)<br> return Anz, Arow, Acol<br><br><br>@njit<br>def define_nonzero_values(<br> rstart:int, <br> rend:int,<br> dx:float,<br> xnodes:int,<br> ynodes:int,<br> Anz:npt.NDArray[np.float64],<br> Arow:npt.NDArray[np.int64],<br> Acol:npt.NDArray[np.int64]<br>) -> int:<br> """ Fill matrix A<br> """<br> icount_nz = 0<br> for row in range(rstart, rend):<br> i, j = index_to_grid(row, xnodes)<br> #A[row, row] = 4.0/dx**2<br> Anz[icount_nz] = 4.0/dx**2<br> Arow[icount_nz] = row<br> Acol[icount_nz] = row<br> icount_nz += 1<br> if i > 0:<br> column = row - xnodes<br> #A[row, column] = -1.0/dx**2<br> Anz[icount_nz] = -1.0/dx**2<br> Arow[icount_nz] = row<br> Acol[icount_nz] = column<br> icount_nz += 1<br> if i < xnodes - 1:<br> column = row + xnodes<br> #A[row, column] = -1.0/dx**2<br> Anz[icount_nz] = -1.0/dx**2<br> Arow[icount_nz] = row<br> Acol[icount_nz] = column<br> icount_nz += 1<br> if j > 0:<br> column = row - 1<br> #A[row, column] = -1.0/dx**2<br> Anz[icount_nz] = -1.0/dx**2<br> Arow[icount_nz] = row<br> Acol[icount_nz] = column<br> icount_nz += 1<br> if j < xnodes - 1:<br> column = row + 1<br> #A[row, column] = -1.0/dx**2<br> Anz[icount_nz] = -1.0/dx**2<br> Arow[icount_nz] = row<br> Acol[icount_nz] = column<br> icount_nz += 1<br> return icount_nz<br><br><br>@njit<br>def clean_nonzero_arrays(<br> icount_nz:int, <br> Anz_ini:npt.NDArray[np.float64], <br> Arow_ini:npt.NDArray[np.float64], <br> Acol_ini:npt.NDArray[np.float64]<br>) -> Tuple[<br> npt.NDArray[np.float64], <br> npt.NDArray[np.float64], <br> npt.NDArray[np.float64]<br> ]:<br> Anz = np.zeros((icount_nz), dtype=np.float64)<br> Arow = np.zeros((icount_nz), dtype=np.int32)<br> Acol = np.zeros((icount_nz), dtype=np.int32)<br> for i in range(icount_nz):<br> Anz[i] = Anz_ini[i]<br> Arow[i] = Arow_ini[i]<br> Acol[i] = Acol_ini[i]<br> return Anz, Arow, Acol<br><br><br>@njit<br>def index_to_grid(r:int, n:int) -> Tuple[int,int]:<br> """Convert a row number into a grid point."""<br> return (r//n, r%n)<br><br><br>if __name__ == "__main__":<br> main()</div><div><br></div></div><div><br></div>
</div><div id="m_-2768164387629248427yahoo_quoted_2889210779">
<div style="font-family:"Helvetica Neue",Helvetica,Arial,sans-serif;font-size:13px;color:rgb(38,40,42)">
<div>
On Friday, August 18, 2023 at 11:24:36 AM CDT, Zhang, Hong <<a href="mailto:hongzhang@anl.gov" target="_blank">hongzhang@anl.gov</a>> wrote:
</div>
<div><br></div>
<div><br></div>
<div><div dir="ltr">You can use this to build a PETSc matrix with the index arrays ai,aj and the value array aa:<br clear="none"> PETSc.Mat().createAIJ(size=(nrows,ncols), csr=(ai,aj,aa))<br clear="none"><br clear="none">Hong (Mr.)<br clear="none"><div id="m_-2768164387629248427yqtfd33436"><br clear="none">> On Aug 17, 2023, at 7:37 AM, Erik Kneller via petsc-users <<a shape="rect" href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a>> wrote:<br clear="none">> <br clear="none">> Hi All,<br clear="none">> <br clear="none">> 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.<br clear="none">> <br clear="none">> 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.<br clear="none">> <br clear="none">> Example Algorithm Building on Lisandro Dalcin's 2D Poisson Example:<br clear="none">> ----------------------------------------------------------------------------------------------<br clear="none">> comm = PETSc.COMM_WORLD<br clear="none">> rank = comm.getRank()<br clear="none">> <br clear="none">> dx = 1.0/(xnodes + 1) # xnodes is the number of nodes in the x and y-directions of the grid<br clear="none">> nnz_max = 5 # max number of non-zero values per row<br clear="none">> <br clear="none">> A = PETSc.Mat()<br clear="none">> A.create(comm=PETSc.COMM_WORLD)<br clear="none">> A.setSizes((xnodes*ynodes, xnodes*ynodes))<br clear="none">> A.setType(PETSc.Mat.Type.AIJ)<br clear="none">> A.setPreallocationNNZ(nnz_max)<br clear="none">> <br clear="none">> rstart, rend = A.getOwnershipRange()<br clear="none">> <br clear="none">> # Here Anz, Arow and Acol are vectors with size equal to the number of non-zero values <br clear="none">> Anz, Arow, Acol = build_nonzero_numpy_arrays_using_numba(rstart, rend, nnz_max, dx, xnodes, ynodes)<br clear="none">> <br clear="none">> A.setValues(Arow, Acol, Anz) # <--- This does not work.<br clear="none">> <br clear="none">> A.assemblyBegin()<br clear="none">> A.assemblyEnd()<br clear="none">> ksp = PETSc.KSP()<br clear="none">> ksp.create(comm=A.getComm())<br clear="none">> ksp.setType(<a href="http://PETSc.KSP.Type.CG" target="_blank">PETSc.KSP.Type.CG</a>)<br clear="none">> ksp.getPC().setType(PETSc.PC.Type.GAMG)<br clear="none">> ksp.setOperators(A)<br clear="none">> ksp.setFromOptions()<br clear="none">> x, b = A.createVecs()<br clear="none">> b.set(1.0)<br clear="none">> <br clear="none">> ksp.solve(b, x)<br clear="none">> <br clear="none">> Regards,<br clear="none">> Erik Kneller, Ph.D.<br clear="none">> <br clear="none"><br clear="none"></div></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>