<html xmlns:v="urn:schemas-microsoft-com:vml" xmlns:o="urn:schemas-microsoft-com:office:office" xmlns:w="urn:schemas-microsoft-com:office:word" xmlns:m="http://schemas.microsoft.com/office/2004/12/omml" xmlns="http://www.w3.org/TR/REC-html40">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=us-ascii">
<meta name="Generator" content="Microsoft Word 14 (filtered medium)">
<style><!--
/* Font Definitions */
@font-face
{font-family:Calibri;
panose-1:2 15 5 2 2 2 4 3 2 4;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0in;
margin-bottom:.0001pt;
font-size:11.0pt;
font-family:"Calibri","sans-serif";}
a:link, span.MsoHyperlink
{mso-style-priority:99;
color:blue;
text-decoration:underline;}
a:visited, span.MsoHyperlinkFollowed
{mso-style-priority:99;
color:purple;
text-decoration:underline;}
span.EmailStyle17
{mso-style-type:personal-compose;
font-family:"Calibri","sans-serif";
color:windowtext;}
.MsoChpDefault
{mso-style-type:export-only;
font-family:"Calibri","sans-serif";}
@page WordSection1
{size:8.5in 11.0in;
margin:1.0in 1.0in 1.0in 1.0in;}
div.WordSection1
{page:WordSection1;}
--></style><!--[if gte mso 9]><xml>
<o:shapedefaults v:ext="edit" spidmax="1026" />
</xml><![endif]--><!--[if gte mso 9]><xml>
<o:shapelayout v:ext="edit">
<o:idmap v:ext="edit" data="1" />
</o:shapelayout></xml><![endif]-->
</head>
<body lang="EN-US" link="blue" vlink="purple">
<div class="WordSection1">
<p class="MsoNormal">Hello,<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">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.
<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">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.<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">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?
<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">Thanks for your help in advance. <o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">Regards,<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">Firat<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">from petsc4py import PETSc<o:p></o:p></p>
<p class="MsoNormal">from mpi4py import MPI<o:p></o:p></p>
<p class="MsoNormal">import numpy as np<o:p></o:p></p>
<p class="MsoNormal">import time<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">comm = MPI.COMM_WORLD<o:p></o:p></p>
<p class="MsoNormal">rank = comm.Get_rank()<o:p></o:p></p>
<p class="MsoNormal">size = comm.Get_size()<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">if rank==0:<o:p></o:p></p>
<p class="MsoNormal"> # proc 0 loads tomo image and does fast calculations to append row, col, val, rh lists<o:p></o:p></p>
<p class="MsoNormal"> # in the real code this vectors will be available on proc 0 no txt files are read<o:p></o:p></p>
<p class="MsoNormal"> row = np.loadtxt('row.out') # indices of non-zero rows <o:p>
</o:p></p>
<p class="MsoNormal"> col = np.loadtxt('col.out') # indices of non-zero columns<o:p></o:p></p>
<p class="MsoNormal"> val = np.loadtxt('vs.out') # values in the sparse matrix<o:p></o:p></p>
<p class="MsoNormal"> rh = np.loadtxt('RHS.out') # right hand side vector<o:p></o:p></p>
<p class="MsoNormal"> n = row.shape[0] #1045699 <o:p></o:p></p>
<p class="MsoNormal"> m = rh.shape[0] #181197 square sparse matrix size<o:p></o:p></p>
<p class="MsoNormal">else: <o:p></o:p></p>
<p class="MsoNormal"> n = None<o:p></o:p></p>
<p class="MsoNormal"> m = None<o:p></o:p></p>
<p class="MsoNormal"> row = None<o:p></o:p></p>
<p class="MsoNormal"> col = None<o:p></o:p></p>
<p class="MsoNormal"> val = None<o:p></o:p></p>
<p class="MsoNormal"> rh = None<o:p></o:p></p>
<p class="MsoNormal"> rh_ind = None<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">m_lcl = comm.bcast(m,root=0) <o:p></o:p></p>
<p class="MsoNormal">n_lcl = comm.bcast(n,root=0)<o:p></o:p></p>
<p class="MsoNormal">neq = n_lcl//size<o:p></o:p></p>
<p class="MsoNormal">meq = m_lcl//size<o:p></o:p></p>
<p class="MsoNormal">nx = np.mod(n_lcl,size)<o:p></o:p></p>
<p class="MsoNormal">mx = np.mod(m_lcl,size)<o:p></o:p></p>
<p class="MsoNormal">row_lcl = np.zeros(neq)<o:p></o:p></p>
<p class="MsoNormal">col_lcl = np.zeros(neq)<o:p></o:p></p>
<p class="MsoNormal">val_lcl = np.zeros(neq)<o:p></o:p></p>
<p class="MsoNormal">rh_lcl = np.zeros(meq)<o:p></o:p></p>
<p class="MsoNormal">a = [neq]*size #send counts for Scatterv<o:p></o:p></p>
<p class="MsoNormal">am = [meq]*size #send counts for Scatterv<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">if nx>0:<o:p></o:p></p>
<p class="MsoNormal"> for i in range(0,nx):<o:p></o:p></p>
<p class="MsoNormal"> if rank==i:<o:p></o:p></p>
<p class="MsoNormal"> row_lcl = np.zeros(neq+1)<o:p></o:p></p>
<p class="MsoNormal"> col_lcl = np.zeros(neq+1)<o:p></o:p></p>
<p class="MsoNormal"> val_lcl = np.zeros(neq+1)<o:p></o:p></p>
<p class="MsoNormal"> a[i] = a[i]+1<o:p></o:p></p>
<p class="MsoNormal">if mx>0:<o:p></o:p></p>
<p class="MsoNormal"> for ii in range(0,mx):<o:p></o:p></p>
<p class="MsoNormal"> if rank==ii:<o:p></o:p></p>
<p class="MsoNormal"> rh_lcl = np.zeros(meq+1)<o:p></o:p></p>
<p class="MsoNormal"> am[ii] = am[ii]+1<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">comm.Scatterv([row,a],row_lcl)<o:p></o:p></p>
<p class="MsoNormal">comm.Scatterv([col,a],col_lcl)<o:p></o:p></p>
<p class="MsoNormal">comm.Scatterv([val,a],val_lcl)<o:p></o:p></p>
<p class="MsoNormal">comm.Scatterv([rh,am],rh_lcl)<o:p></o:p></p>
<p class="MsoNormal">comm.Barrier()<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">A = PETSc.Mat()<o:p></o:p></p>
<p class="MsoNormal">A.create()<o:p></o:p></p>
<p class="MsoNormal">A.setSizes([m_lcl,m_lcl])<o:p></o:p></p>
<p class="MsoNormal">A.setType('aij')<o:p></o:p></p>
<p class="MsoNormal">A.setUp()<o:p></o:p></p>
<p class="MsoNormal">lr = row_lcl.shape[0]<o:p></o:p></p>
<p class="MsoNormal">for i in range(0,lr):<o:p></o:p></p>
<p class="MsoNormal"> A[row_lcl[i],col_lcl[i]] = val_lcl[i]<o:p></o:p></p>
<p class="MsoNormal">A.assemblyBegin()<o:p></o:p></p>
<p class="MsoNormal">A.assemblyEnd()<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">if size>1: # to get the range for scattered vectors<o:p></o:p></p>
<p class="MsoNormal"> ami = [0]<o:p></o:p></p>
<p class="MsoNormal"> ami = np.array([0]+am).cumsum()<o:p></o:p></p>
<p class="MsoNormal"> for kk in range(0,size):<o:p></o:p></p>
<p class="MsoNormal"> if rank==kk:<o:p></o:p></p>
<p class="MsoNormal"> Is = ami[kk]<o:p></o:p></p>
<p class="MsoNormal"> Ie = ami[kk+1]<o:p></o:p></p>
<p class="MsoNormal">else:<o:p></o:p></p>
<p class="MsoNormal"> Is=0; Ie=m_lcl<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">b= PETSc.Vec()<o:p></o:p></p>
<p class="MsoNormal">b.create()<o:p></o:p></p>
<p class="MsoNormal">b.setSizes(m_lcl)<o:p></o:p></p>
<p class="MsoNormal">b.setFromOptions() <o:p></o:p></p>
<p class="MsoNormal">b.setUp()<o:p></o:p></p>
<p class="MsoNormal">b.setValues(list(range(Is,Ie)),rh_lcl)<o:p></o:p></p>
<p class="MsoNormal">b.assemblyBegin()<o:p></o:p></p>
<p class="MsoNormal">b.assemblyEnd()<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"># solution vector<o:p></o:p></p>
<p class="MsoNormal">x = b.duplicate() <o:p></o:p></p>
<p class="MsoNormal">x.assemblyBegin()<o:p></o:p></p>
<p class="MsoNormal">x.assemblyEnd()<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"># create linear solver<o:p></o:p></p>
<p class="MsoNormal">ksp = PETSc.KSP()<o:p></o:p></p>
<p class="MsoNormal">ksp.create()<o:p></o:p></p>
<p class="MsoNormal">ksp.setOperators(A)<o:p></o:p></p>
<p class="MsoNormal">ksp.setType('cg')<o:p></o:p></p>
<p class="MsoNormal">#ksp.getPC().setType('icc') # only sequential<o:p></o:p></p>
<p class="MsoNormal">ksp.getPC().setType('jacobi')<o:p></o:p></p>
<p class="MsoNormal">print('solving with:', ksp.getType())<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">#solve<o:p></o:p></p>
<p class="MsoNormal">st=time.time()<o:p></o:p></p>
<p class="MsoNormal">ksp.solve(b,x)<o:p></o:p></p>
<p class="MsoNormal">et=time.time()<o:p></o:p></p>
<p class="MsoNormal">print(et-st)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">if size>1:<o:p></o:p></p>
<p class="MsoNormal"> #gather<o:p></o:p></p>
<p class="MsoNormal"> if rank==0:<o:p></o:p></p>
<p class="MsoNormal"> xGthr = np.zeros(m)<o:p></o:p></p>
<p class="MsoNormal"> else:<o:p></o:p></p>
<p class="MsoNormal"> xGthr = None <o:p></o:p></p>
<p class="MsoNormal"> comm.Gatherv(x,[xGthr,am])<o:p></o:p></p>
</div>
</body>
</html>