[petsc-users] petsc4py sparse matrix construction time

[Matthew Knepley]

On Mon, Oct 30, 2017 at 8:06 PM, Cetinbas, Cankur Firat <ccetinbas at anl.gov>
wrote:

> Hello,
>
>
>
> I am a beginner both in PETSc and mpi4py. I have been working on
> parallelizing our water transport code (where we solve linear system of
> equations) and I started with the toy code below.
>
>
>
> The toy code reads right hand size (rh), row, column, value vectors to
> construct sparse coefficient matrix and then scatters them to construct
> parallel PETSc coefficient matrix and right hand side vector.
>
>
>
> The sparse matrix generation time is extremely high in comparison to
> sps.csr_matrix((val, (row, col)), shape=(n,n)) in python. For instance
> python generates 181197x181197 sparse matrix in 0.06 seconds and this code
> with 32 cores:1.19s, 16 cores:6.98s and 8 cores:29.5 s.  I was wondering if
> I am making a mistake in generating sparse matrix? Is there a more
> efficient way?
>

It looks like you do not preallocate the matrix. There is a chapter on this
in the manual.

   Matt


> Thanks for your help in advance.
>
>
>
> Regards,
>
>
>
> Firat
>
>
>
> from petsc4py import PETSc
>
> from mpi4py import MPI
>
> import numpy as np
>
> import time
>
>
>
> comm = MPI.COMM_WORLD
>
> rank = comm.Get_rank()
>
> size = comm.Get_size()
>
>
>
> if rank==0:
>
>     # proc 0 loads tomo image and does fast calculations to append row,
> col, val, rh lists
>
>     # in the real code this vectors will be available on proc 0 no txt
> files are read
>
>     row = np.loadtxt('row.out')  # indices of non-zero rows
>
>     col = np.loadtxt('col.out') # indices of non-zero columns
>
>     val = np.loadtxt('vs.out')  # values in the sparse matrix
>
>     rh = np.loadtxt('RHS.out') # right hand side vector
>
>     n = row.shape[0]  #1045699
>
>     m = rh.shape[0] #181197 square sparse matrix size
>
> else:
>
>     n = None
>
>     m = None
>
>     row = None
>
>     col = None
>
>     val = None
>
>     rh = None
>
>     rh_ind = None
>
>
>
> m_lcl = comm.bcast(m,root=0)
>
> n_lcl = comm.bcast(n,root=0)
>
> neq = n_lcl//size
>
> meq = m_lcl//size
>
> nx = np.mod(n_lcl,size)
>
> mx = np.mod(m_lcl,size)
>
> row_lcl = np.zeros(neq)
>
> col_lcl = np.zeros(neq)
>
> val_lcl = np.zeros(neq)
>
> rh_lcl = np.zeros(meq)
>
> a = [neq]*size #send counts for Scatterv
>
> am = [meq]*size #send counts for Scatterv
>
>
>
> if nx>0:
>
>     for i in range(0,nx):
>
>         if rank==i:
>
>             row_lcl = np.zeros(neq+1)
>
>             col_lcl = np.zeros(neq+1)
>
>             val_lcl = np.zeros(neq+1)
>
>         a[i] = a[i]+1
>
> if mx>0:
>
>     for ii in range(0,mx):
>
>         if rank==ii:
>
>             rh_lcl = np.zeros(meq+1)
>
>         am[ii] = am[ii]+1
>
>
>
> comm.Scatterv([row,a],row_lcl)
>
> comm.Scatterv([col,a],col_lcl)
>
> comm.Scatterv([val,a],val_lcl)
>
> comm.Scatterv([rh,am],rh_lcl)
>
> comm.Barrier()
>
>
>
> A = PETSc.Mat()
>
> A.create()
>
> A.setSizes([m_lcl,m_lcl])
>
> A.setType('aij')
>
> A.setUp()
>
> lr = row_lcl.shape[0]
>
> for i in range(0,lr):
>
>     A[row_lcl[i],col_lcl[i]] = val_lcl[i]
>
> A.assemblyBegin()
>
> A.assemblyEnd()
>
>
>
> if size>1: # to get the range for scattered vectors
>
>     ami = [0]
>
>     ami = np.array([0]+am).cumsum()
>
>     for kk in range(0,size):
>
>         if rank==kk:
>
>             Is = ami[kk]
>
>             Ie = ami[kk+1]
>
> else:
>
>     Is=0; Ie=m_lcl
>
>
>
> b= PETSc.Vec()
>
> b.create()
>
> b.setSizes(m_lcl)
>
> b.setFromOptions()
>
> b.setUp()
>
> b.setValues(list(range(Is,Ie)),rh_lcl)
>
> b.assemblyBegin()
>
> b.assemblyEnd()
>
>
>
> # solution vector
>
> x = b.duplicate()
>
> x.assemblyBegin()
>
> x.assemblyEnd()
>
>
>
> # create linear solver
>
> ksp = PETSc.KSP()
>
> ksp.create()
>
> ksp.setOperators(A)
>
> ksp.setType('cg')
>
> #ksp.getPC().setType('icc') # only sequential
>
> ksp.getPC().setType('jacobi')
>
> print('solving with:', ksp.getType())
>
>
>
> #solve
>
> st=time.time()
>
> ksp.solve(b,x)
>
> et=time.time()
>
> print(et-st)
>
>
>
> if size>1:
>
>     #gather
>
>     if rank==0:
>
>         xGthr = np.zeros(m)
>
>     else:
>
>         xGthr = None
>
>     comm.Gatherv(x,[xGthr,am])
>



-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ <http://www.caam.rice.edu/~mk51/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20171030/cfa114c9/attachment-0001.html>