[petsc-users] petsc4py sparse matrix construction time

[Cetinbas, Cankur Firat]

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?

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])
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20171031/5b5bb841/attachment.html>