<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:x="urn:schemas-microsoft-com:office:excel" 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=utf-8">
<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;}
@font-face
{font-family:Tahoma;
panose-1:2 11 6 4 3 5 4 4 2 4;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0in;
margin-bottom:.0001pt;
font-size:12.0pt;
font-family:"Times New Roman","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;}
p.MsoAcetate, li.MsoAcetate, div.MsoAcetate
{mso-style-priority:99;
mso-style-link:"Balloon Text Char";
margin:0in;
margin-bottom:.0001pt;
font-size:8.0pt;
font-family:"Tahoma","sans-serif";}
p.MsoListParagraph, li.MsoListParagraph, div.MsoListParagraph
{mso-style-priority:34;
margin-top:0in;
margin-right:0in;
margin-bottom:0in;
margin-left:.5in;
margin-bottom:.0001pt;
font-size:12.0pt;
font-family:"Times New Roman","serif";}
span.hoenzb
{mso-style-name:hoenzb;}
span.EmailStyle18
{mso-style-type:personal-reply;
font-family:"Calibri","sans-serif";
color:#1F497D;}
span.BalloonTextChar
{mso-style-name:"Balloon Text Char";
mso-style-priority:99;
mso-style-link:"Balloon Text";
font-family:"Tahoma","sans-serif";}
.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;}
/* List Definitions */
@list l0
{mso-list-id:1897623605;
mso-list-type:hybrid;
mso-list-template-ids:1420994276 -1875214254 67698713 67698715 67698703 67698713 67698715 67698703 67698713 67698715;}
@list l0:level1
{mso-level-text:"%1\)";
mso-level-tab-stop:none;
mso-level-number-position:left;
margin-left:24.0pt;
text-indent:-.25in;}
@list l0:level2
{mso-level-number-format:alpha-lower;
mso-level-tab-stop:none;
mso-level-number-position:left;
margin-left:60.0pt;
text-indent:-.25in;}
@list l0:level3
{mso-level-number-format:roman-lower;
mso-level-tab-stop:none;
mso-level-number-position:right;
margin-left:96.0pt;
text-indent:-9.0pt;}
@list l0:level4
{mso-level-tab-stop:none;
mso-level-number-position:left;
margin-left:132.0pt;
text-indent:-.25in;}
@list l0:level5
{mso-level-number-format:alpha-lower;
mso-level-tab-stop:none;
mso-level-number-position:left;
margin-left:168.0pt;
text-indent:-.25in;}
@list l0:level6
{mso-level-number-format:roman-lower;
mso-level-tab-stop:none;
mso-level-number-position:right;
margin-left:204.0pt;
text-indent:-9.0pt;}
@list l0:level7
{mso-level-tab-stop:none;
mso-level-number-position:left;
margin-left:240.0pt;
text-indent:-.25in;}
@list l0:level8
{mso-level-number-format:alpha-lower;
mso-level-tab-stop:none;
mso-level-number-position:left;
margin-left:276.0pt;
text-indent:-.25in;}
@list l0:level9
{mso-level-number-format:roman-lower;
mso-level-tab-stop:none;
mso-level-number-position:right;
margin-left:312.0pt;
text-indent:-9.0pt;}
ol
{margin-bottom:0in;}
ul
{margin-bottom:0in;}
--></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"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D">It took me a while to respond because I couldn’t get converged results with the iterative solvers (they didn’t give any errors but the results didn’t match
with Matlab mldivide and Python spsolve). So we spent some time trying to reconfigure petsc and petsc4py with MUMPS.
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D">I don’t know if locally called mat.createAIJ or setValues come with any overhead but I timed each section of the code locally.
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D"><o:p> </o:p></span></p>
<table class="MsoNormalTable" border="0" cellspacing="0" cellpadding="0" width="494" style="width:370.15pt;margin-left:-.75pt;border-collapse:collapse">
<tbody>
<tr style="height:15.0pt">
<td width="127" nowrap="" valign="bottom" style="width:95.25pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
</td>
<td width="91" nowrap="" valign="bottom" style="width:68.15pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal"><b><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">Single core<o:p></o:p></span></b></p>
</td>
<td width="64" nowrap="" valign="bottom" style="width:48.0pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal"><b><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">2 cores<o:p></o:p></span></b></p>
</td>
<td width="64" nowrap="" valign="bottom" style="width:48.0pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal"><b><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">4 cores
<o:p></o:p></span></b></p>
</td>
<td width="65" nowrap="" valign="bottom" style="width:48.5pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal"><b><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">8 cores
<o:p></o:p></span></b></p>
</td>
<td width="83" nowrap="" valign="bottom" style="width:62.25pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal"><b><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">16 cores<o:p></o:p></span></b></p>
</td>
</tr>
<tr style="height:15.0pt">
<td width="127" nowrap="" valign="bottom" style="width:95.25pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal"><b><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">bcast<o:p></o:p></span></b></p>
</td>
<td width="91" nowrap="" valign="bottom" style="width:68.15pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
</td>
<td width="64" nowrap="" valign="bottom" style="width:48.0pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">0.064<o:p></o:p></span></p>
</td>
<td width="64" nowrap="" valign="bottom" style="width:48.0pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">0.064<o:p></o:p></span></p>
</td>
<td width="65" nowrap="" valign="bottom" style="width:48.5pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">0.071<o:p></o:p></span></p>
</td>
<td width="83" nowrap="" valign="bottom" style="width:62.25pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">0.069<o:p></o:p></span></p>
</td>
</tr>
<tr style="height:15.0pt">
<td width="127" nowrap="" valign="bottom" style="width:95.25pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal"><b><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">local vector init<o:p></o:p></span></b></p>
</td>
<td width="91" nowrap="" valign="bottom" style="width:68.15pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
</td>
<td width="64" nowrap="" valign="bottom" style="width:48.0pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">0.004<o:p></o:p></span></p>
</td>
<td width="64" nowrap="" valign="bottom" style="width:48.0pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">0.0028<o:p></o:p></span></p>
</td>
<td width="65" nowrap="" valign="bottom" style="width:48.5pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">0.0015<o:p></o:p></span></p>
</td>
<td width="83" nowrap="" valign="bottom" style="width:62.25pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">0.00087<o:p></o:p></span></p>
</td>
</tr>
<tr style="height:15.0pt">
<td width="127" nowrap="" valign="bottom" style="width:95.25pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal"><b><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">scatterv<o:p></o:p></span></b></p>
</td>
<td width="91" nowrap="" valign="bottom" style="width:68.15pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
</td>
<td width="64" nowrap="" valign="bottom" style="width:48.0pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">0.0076<o:p></o:p></span></p>
</td>
<td width="64" nowrap="" valign="bottom" style="width:48.0pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">0.0056<o:p></o:p></span></p>
</td>
<td width="65" nowrap="" valign="bottom" style="width:48.5pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">0.0079<o:p></o:p></span></p>
</td>
<td width="83" nowrap="" valign="bottom" style="width:62.25pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">0.0099<o:p></o:p></span></p>
</td>
</tr>
<tr style="height:15.0pt">
<td width="127" nowrap="" valign="bottom" style="width:95.25pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal"><b><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">convert to int 32<o:p></o:p></span></b></p>
</td>
<td width="91" nowrap="" valign="bottom" style="width:68.15pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
</td>
<td width="64" nowrap="" valign="bottom" style="width:48.0pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">0.0024<o:p></o:p></span></p>
</td>
<td width="64" nowrap="" valign="bottom" style="width:48.0pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">0.0015<o:p></o:p></span></p>
</td>
<td width="65" nowrap="" valign="bottom" style="width:48.5pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">0.00094<o:p></o:p></span></p>
</td>
<td width="83" nowrap="" valign="bottom" style="width:62.25pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">0.00072<o:p></o:p></span></p>
</td>
</tr>
<tr style="height:15.0pt">
<td width="127" nowrap="" valign="bottom" style="width:95.25pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal"><b><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">calculate indptr<o:p></o:p></span></b></p>
</td>
<td width="91" nowrap="" valign="bottom" style="width:68.15pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
</td>
<td width="64" nowrap="" valign="bottom" style="width:48.0pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">0.0028<o:p></o:p></span></p>
</td>
<td width="64" nowrap="" valign="bottom" style="width:48.0pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">0.0018<o:p></o:p></span></p>
</td>
<td width="65" nowrap="" valign="bottom" style="width:48.5pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">0.0011<o:p></o:p></span></p>
</td>
<td width="83" nowrap="" valign="bottom" style="width:62.25pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">0.0007<o:p></o:p></span></p>
</td>
</tr>
<tr style="height:15.0pt">
<td width="127" nowrap="" valign="bottom" style="width:95.25pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal"><b><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">create AIJ<o:p></o:p></span></b></p>
</td>
<td width="91" nowrap="" valign="bottom" style="width:68.15pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">0.044<o:p></o:p></span></p>
</td>
<td width="64" nowrap="" valign="bottom" style="width:48.0pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">0.0037<o:p></o:p></span></p>
</td>
<td width="64" nowrap="" valign="bottom" style="width:48.0pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">0.022<o:p></o:p></span></p>
</td>
<td width="65" nowrap="" valign="bottom" style="width:48.5pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">0.017<o:p></o:p></span></p>
</td>
<td width="83" nowrap="" valign="bottom" style="width:62.25pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">0.013<o:p></o:p></span></p>
</td>
</tr>
<tr style="height:15.0pt">
<td width="127" nowrap="" valign="bottom" style="width:95.25pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal"><b><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">RHS vector<o:p></o:p></span></b></p>
</td>
<td width="91" nowrap="" valign="bottom" style="width:68.15pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
</td>
<td width="64" nowrap="" valign="bottom" style="width:48.0pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">0.0016<o:p></o:p></span></p>
</td>
<td width="64" nowrap="" valign="bottom" style="width:48.0pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">0.0017<o:p></o:p></span></p>
</td>
<td width="65" nowrap="" valign="bottom" style="width:48.5pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">0.0019<o:p></o:p></span></p>
</td>
<td width="83" nowrap="" valign="bottom" style="width:62.25pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">0.0023<o:p></o:p></span></p>
</td>
</tr>
<tr style="height:15.0pt">
<td width="127" nowrap="" valign="bottom" style="width:95.25pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal"><b><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">solve (MUMPS)<o:p></o:p></span></b></p>
</td>
<td width="91" nowrap="" valign="bottom" style="width:68.15pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">10.94<o:p></o:p></span></p>
</td>
<td width="64" nowrap="" valign="bottom" style="width:48.0pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">7.23<o:p></o:p></span></p>
</td>
<td width="64" nowrap="" valign="bottom" style="width:48.0pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">5.67<o:p></o:p></span></p>
</td>
<td width="65" nowrap="" valign="bottom" style="width:48.5pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">5.78<o:p></o:p></span></p>
</td>
<td width="83" nowrap="" valign="bottom" style="width:62.25pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">5.13<o:p></o:p></span></p>
</td>
</tr>
<tr style="height:15.0pt">
<td width="127" nowrap="" valign="bottom" style="width:95.25pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal"><b><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">total time (s)<o:p></o:p></span></b></p>
</td>
<td width="91" nowrap="" valign="bottom" style="width:68.15pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">10.984<o:p></o:p></span></p>
</td>
<td width="64" nowrap="" valign="bottom" style="width:48.0pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">7.3161<o:p></o:p></span></p>
</td>
<td width="64" nowrap="" valign="bottom" style="width:48.0pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">5.7694<o:p></o:p></span></p>
</td>
<td width="65" nowrap="" valign="bottom" style="width:48.5pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">5.88134<o:p></o:p></span></p>
</td>
<td width="83" nowrap="" valign="bottom" style="width:62.25pt;padding:0in 5.4pt 0in 5.4pt;height:15.0pt">
<p class="MsoNormal" align="right" style="text-align:right"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:black">5.22649<o:p></o:p></span></p>
</td>
</tr>
</tbody>
</table>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D">MUMPS solver is the most time consuming part so I need to find a way to make iterative solvers converge. The single core MUMPS solver takes almost 9 times longer
with respect to mldivide in Matlab. I might be doing something wrong setting up the solver but the results of MUMPS match with Matlab.
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D"><o:p> </o:p></span></p>
<p class="MsoNormal" style="text-autospace:none"><span style="font-size:11.0pt;font-family:"Courier New""> ksp = PETSc.KSP()<o:p></o:p></span></p>
<p class="MsoNormal" style="text-autospace:none"><span style="font-size:11.0pt;font-family:"Courier New""> ksp.create()<o:p></o:p></span></p>
<p class="MsoNormal" style="text-autospace:none"><span style="font-size:11.0pt;font-family:"Courier New""> ksp.setOperators(pA)<o:p></o:p></span></p>
<p class="MsoNormal" style="text-autospace:none"><span style="font-size:11.0pt;font-family:"Courier New""> ksp.setType('preonly')<o:p></o:p></span></p>
<p class="MsoNormal" style="text-autospace:none"><span style="font-size:11.0pt;font-family:"Courier New""> pc = ksp.getPC()<o:p></o:p></span></p>
<p class="MsoNormal" style="text-autospace:none"><span style="font-size:11.0pt;font-family:"Courier New""> pc.setType('lu')<o:p></o:p></span></p>
<p class="MsoNormal" style="text-autospace:none"><span style="font-size:11.0pt;font-family:"Courier New""> pc.setFactorSolverPackage('mumps')<o:p></o:p></span></p>
<p class="MsoNormal" style="text-autospace:none"><span style="font-size:11.0pt;font-family:"Courier New""> ksp.solve(rhv,x)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D">In addition I don’t understand why broadcast takes so much time w.r.t. scatter. I only pass a small array at the size of number of processes.<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D">Regards,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D">Firat<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D"><o:p> </o:p></span></p>
<p class="MsoNormal"><b><span style="font-size:10.0pt;font-family:"Tahoma","sans-serif"">From:</span></b><span style="font-size:10.0pt;font-family:"Tahoma","sans-serif""> Matthew Knepley [mailto:knepley@gmail.com]
<br>
<b>Sent:</b> Wednesday, November 01, 2017 3:35 AM<br>
<b>To:</b> Cetinbas, Cankur Firat<br>
<b>Cc:</b> Smith, Barry F.; petsc-users@mcs.anl.gov; Ahluwalia, Rajesh K.<br>
<b>Subject:</b> Re: [petsc-users] petsc4py sparse matrix construction time<o:p></o:p></span></p>
<p class="MsoNormal"><o:p> </o:p></p>
<div>
<div>
<div>
<p class="MsoNormal">On Tue, Oct 31, 2017 at 10:00 PM, Cetinbas, Cankur Firat <<a href="mailto:ccetinbas@anl.gov" target="_blank">ccetinbas@anl.gov</a>> wrote:<o:p></o:p></p>
<p class="MsoNormal">Hi,<br>
<br>
Thanks a lot. Based on both of your suggestions I modified the code using Mat.createAIJ() and csr option. The computation time decreased significantly after using this method. Still if there is a better option please let me know after seeing the modified code
below.<br>
<br>
At first trial with 1000x1000 matrix with 96019 non-zeros in the matrix, the computation time did not scale with the number of cores : Single core python @ 0.0035s, single core petsc @ 0.0024s, 2 cores petsc @ 0.0036s, 4 cores petsc @ 0.0032, 8 cores petsc
@ 0.0030s.<br>
<br>
Then I tried with larger matrix 181797x181797 with more non-zeros and I got the following results: Single core python @ 0.021, single core petsc @ 0.031, 2 cores petsc @ 0.024s, 4 cores petsc @ 0.014, 8 cores petsc @ 0.009s, 16 cores petsc @ 0.0087s.<br>
<br>
I think the optimum number of nodes is highly dependent on matrix size and the number of non-zeros. In the real code matrix size (and so the number of non-zero elements) will grow at every iteration starting with very small matrices growing to very big ones.
Is it possible to set the number process from the code dynamically?<o:p></o:p></p>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">I am not sure how you are interpreting these measurements. Normally, I would say<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoListParagraph" style="margin-left:24.0pt;text-indent:-.25in;mso-list:l0 level1 lfo1">
<![if !supportLists]><span style="color:#1F497D"><span style="mso-list:Ignore">1)<span style="font:7.0pt "Times New Roman"">
</span></span></span><![endif]>Everything timed below is "parallel overhead". This is not intended to scale with P, instead it will look like a constant, as you observe<span style="color:#1F497D"><o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal"> 2) The time to compute the matrix entires should far outstrip the time below to figure the nonzero structure. Is this true?<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal"> 3) Solve time is often larger than matrix calculation. Is it?<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">Thus, we deciding on parallelism, you need to look at the largest costs, and how they scale with P.<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal"> Thanks,<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal"> Matt<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<blockquote style="border:none;border-left:solid #CCCCCC 1.0pt;padding:0in 0in 0in 6.0pt;margin-left:4.8pt;margin-right:0in">
<p class="MsoNormal" style="margin-bottom:12.0pt">Another question is about the data types; mpi4py only let me transfer float type data, and petsc4py only lets me use int32 type indices. Besides keep converting the data, is there any solution for this?<br>
<br>
The modified code for matrix creation part:<br>
<br>
comm = MPI.COMM_WORLD<br>
rank = comm.Get_rank()<br>
size = comm.Get_size()<br>
<br>
if rank==0:<br>
row = np.loadtxt('row1000.out').astype(dtype='int32')<br>
col = np.loadtxt('col1000.out').astype(dtype='int32')<br>
val = np.loadtxt('val1000.out').astype(dtype='int32')<br>
m = 1000 # 1000 x 1000 matrix<br>
if size>1:<br>
rbc = np.bincount(row)*1.0<br>
ieq = int(np.floor(m/size))<br>
a = [ieq]*size<br>
ix = int(np.mod(m,size))<br>
if ix>0:<br>
for i in range(0,ix):<br>
a[i]= a[i]+1<br>
a = np.array([0]+a).cumsum()<br>
b = np.zeros(a.shape[0]-1)<br>
for i in range(0,a.shape[0]-1):<br>
b[i]=rbc[a[i]:a[i+1]].sum() # b is the send counts for Scatterv<br>
row = row.astype(dtype=float)<br>
col = col.astype(dtype=float)<br>
val = val.astype(dtype=float)<br>
else:<br>
row=None<br>
col=None<br>
val=None<br>
indpt=None<br>
b=None<br>
m=None<br>
<br>
if size>1:<br>
ml = comm.bcast(m,root=0)<br>
bl = comm.bcast(b,root=0)<br>
row_lcl = np.zeros(bl[rank])<br>
col_lcl = row_lcl.copy()<br>
val_lcl = row_lcl.copy()<br>
<br>
comm.Scatterv([row,b],row_lcl)<br>
comm.Scatterv([col,b],col_lcl)<br>
comm.Scatterv([val,b],val_lcl)<br>
comm.Barrier()<br>
<br>
row_lcl = row_lcl.astype(dtype='int32')<br>
col_lcl = col_lcl.astype(dtype='int32')<br>
val_lcl = val_lcl.astype(dtype='int32')<br>
<br>
indptr = np.bincount(row_lcl)<br>
indptr = indptr[indptr>0]<br>
indptr = np.insert(indptr,0,0).cumsum()<br>
indptr = indptr.astype(dtype='int32')<br>
comm.Barrier()<br>
<br>
pA = PETSc.Mat().createAIJ(size=(ml,ml),csr=(indptr, col_lcl, val_lcl)) # Matrix generation<br>
<br>
else:<br>
indptr = np.bincount(row)<br>
indptr = np.insert(indptr,0,0).cumsum()<br>
indptr = indptr.astype(dtype='int32')<br>
st=time.time()<br>
pA = PETSc.Mat().createAIJ(size=(m,m),csr=(indptr, col, val))<br>
print('dt:',time.time()-st)<br>
<br>
<br>
Regards,<br>
<br>
Firat<br>
<br>
<br>
-----Original Message-----<br>
From: Smith, Barry F.<br>
Sent: Tuesday, October 31, 2017 10:18 AM<br>
To: Matthew Knepley<br>
Cc: Cetinbas, Cankur Firat; <a href="mailto:petsc-users@mcs.anl.gov">petsc-users@mcs.anl.gov</a>; Ahluwalia, Rajesh K.<br>
Subject: Re: [petsc-users] petsc4py sparse matrix construction time<br>
<br>
<br>
You also need to make sure that most matrix entries are generated on the process that they will belong on.<br>
<br>
Barry<br>
<br>
> On Oct 30, 2017, at 8:01 PM, Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>> wrote:<br>
><br>
> On Mon, Oct 30, 2017 at 8:06 PM, Cetinbas, Cankur Firat <<a href="mailto:ccetinbas@anl.gov">ccetinbas@anl.gov</a>> wrote:<br>
> Hello,<br>
><br>
><br>
><br>
> 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.<br>
><br>
><br>
><br>
> 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.<br>
><br>
><br>
><br>
> 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?<br>
><br>
><br>
> It looks like you do not preallocate the matrix. There is a chapter on this in the manual.<br>
><br>
> Matt<br>
><br>
> Thanks for your help in advance.<br>
><br>
><br>
><br>
> Regards,<br>
><br>
><br>
><br>
> Firat<br>
><br>
><br>
><br>
> from petsc4py import PETSc<br>
><br>
> from mpi4py import MPI<br>
><br>
> import numpy as np<br>
><br>
> import time<br>
><br>
><br>
><br>
> comm = MPI.COMM_WORLD<br>
><br>
> rank = comm.Get_rank()<br>
><br>
> size = comm.Get_size()<br>
><br>
><br>
><br>
> if rank==0:<br>
><br>
> # proc 0 loads tomo image and does fast calculations to append row, col, val, rh lists<br>
><br>
> # in the real code this vectors will be available on proc 0 no txt files are read<br>
><br>
> row = np.loadtxt('row.out') # indices of non-zero rows<br>
><br>
> col = np.loadtxt('col.out') # indices of non-zero columns<br>
><br>
> val = np.loadtxt('vs.out') # values in the sparse matrix<br>
><br>
> rh = np.loadtxt('RHS.out') # right hand side vector<br>
><br>
> n = row.shape[0] #1045699<br>
><br>
> m = rh.shape[0] #181197 square sparse matrix size<br>
><br>
> else:<br>
><br>
> n = None<br>
><br>
> m = None<br>
><br>
> row = None<br>
><br>
> col = None<br>
><br>
> val = None<br>
><br>
> rh = None<br>
><br>
> rh_ind = None<br>
><br>
><br>
><br>
> m_lcl = comm.bcast(m,root=0)<br>
><br>
> n_lcl = comm.bcast(n,root=0)<br>
><br>
> neq = n_lcl//size<br>
><br>
> meq = m_lcl//size<br>
><br>
> nx = np.mod(n_lcl,size)<br>
><br>
> mx = np.mod(m_lcl,size)<br>
><br>
> row_lcl = np.zeros(neq)<br>
><br>
> col_lcl = np.zeros(neq)<br>
><br>
> val_lcl = np.zeros(neq)<br>
><br>
> rh_lcl = np.zeros(meq)<br>
><br>
> a = [neq]*size #send counts for Scatterv<br>
><br>
> am = [meq]*size #send counts for Scatterv<br>
><br>
><br>
><br>
> if nx>0:<br>
><br>
> for i in range(0,nx):<br>
><br>
> if rank==i:<br>
><br>
> row_lcl = np.zeros(neq+1)<br>
><br>
> col_lcl = np.zeros(neq+1)<br>
><br>
> val_lcl = np.zeros(neq+1)<br>
><br>
> a[i] = a[i]+1<br>
><br>
> if mx>0:<br>
><br>
> for ii in range(0,mx):<br>
><br>
> if rank==ii:<br>
><br>
> rh_lcl = np.zeros(meq+1)<br>
><br>
> am[ii] = am[ii]+1<br>
><br>
><br>
><br>
> comm.Scatterv([row,a],row_lcl)<br>
><br>
> comm.Scatterv([col,a],col_lcl)<br>
><br>
> comm.Scatterv([val,a],val_lcl)<br>
><br>
> comm.Scatterv([rh,am],rh_lcl)<br>
><br>
> comm.Barrier()<br>
><br>
><br>
><br>
> A = PETSc.Mat()<br>
><br>
> A.create()<br>
><br>
> A.setSizes([m_lcl,m_lcl])<br>
><br>
> A.setType('aij')<br>
><br>
> A.setUp()<br>
><br>
> lr = row_lcl.shape[0]<br>
><br>
> for i in range(0,lr):<br>
><br>
> A[row_lcl[i],col_lcl[i]] = val_lcl[i]<br>
><br>
> A.assemblyBegin()<br>
><br>
> A.assemblyEnd()<br>
><br>
><br>
><br>
> if size>1: # to get the range for scattered vectors<br>
><br>
> ami = [0]<br>
><br>
> ami = np.array([0]+am).cumsum()<br>
><br>
> for kk in range(0,size):<br>
><br>
> if rank==kk:<br>
><br>
> Is = ami[kk]<br>
><br>
> Ie = ami[kk+1]<br>
><br>
> else:<br>
><br>
> Is=0; Ie=m_lcl<br>
><br>
><br>
><br>
> b= PETSc.Vec()<br>
><br>
> b.create()<br>
><br>
> b.setSizes(m_lcl)<br>
><br>
> b.setFromOptions()<br>
><br>
> b.setUp()<br>
><br>
> b.setValues(list(range(Is,Ie)),rh_lcl)<br>
><br>
> b.assemblyBegin()<br>
><br>
> b.assemblyEnd()<br>
><br>
><br>
><br>
> # solution vector<br>
><br>
> x = b.duplicate()<br>
><br>
> x.assemblyBegin()<br>
><br>
> x.assemblyEnd()<br>
><br>
><br>
><br>
> # create linear solver<br>
><br>
> ksp = PETSc.KSP()<br>
><br>
> ksp.create()<br>
><br>
> ksp.setOperators(A)<br>
><br>
> ksp.setType('cg')<br>
><br>
> #ksp.getPC().setType('icc') # only sequential<br>
><br>
> ksp.getPC().setType('jacobi')<br>
><br>
> print('solving with:', ksp.getType())<br>
><br>
><br>
><br>
> #solve<br>
><br>
> st=time.time()<br>
><br>
> ksp.solve(b,x)<br>
><br>
> et=time.time()<br>
><br>
> print(et-st)<br>
><br>
><br>
><br>
> if size>1:<br>
><br>
> #gather<br>
><br>
> if rank==0:<br>
><br>
> xGthr = np.zeros(m)<br>
><br>
> else:<br>
><br>
> xGthr = None<br>
><br>
> comm.Gatherv(x,[xGthr,am])<br>
><br>
><br>
><br>
<span class="hoenzb"><span style="color:#888888">></span></span><span style="color:#888888"><br>
<span class="hoenzb">> --</span><br>
<span class="hoenzb">> What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.</span><br>
<span class="hoenzb">> -- Norbert Wiener</span><br>
<span class="hoenzb">></span><br>
<span class="hoenzb">> <a href="https://www.cse.buffalo.edu/~knepley/" target="_blank">
https://www.cse.buffalo.edu/~knepley/</a></span></span><o:p></o:p></p>
</blockquote>
</div>
<p class="MsoNormal"><br>
<br clear="all">
<o:p></o:p></p>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<p class="MsoNormal">-- <o:p></o:p></p>
<div>
<div>
<div>
<div>
<div>
<p class="MsoNormal">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<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal"><a href="http://www.caam.rice.edu/~mk51/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><o:p></o:p></p>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</body>
</html>