<div dir="ltr"><div dir="ltr">On Thu, Mar 11, 2021 at 3:27 AM Mathieu Dutour <<a href="mailto:mathieu.dutour@gmail.com">mathieu.dutour@gmail.com</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 dir="ltr">Dear all,<div><br></div><div>I would like to work with a special kind of linear system that ought to be</div><div>very common but I am not sure that it is possible in PETSC.</div><div><br></div><div>What we have is an unstructured grid with say 3.10^5 nodes in it.</div><div>At each node, we have a number of frequency/direction and together</div><div>this makes about 1000 values at the node. So, in total the linear system</div><div>has say 3.10^8 values.</div><div><br></div><div>We managed to implement this system with Petsc but the performance</div><div>was unsatisfactory. We think that Petsc is not exploiting the special</div><div>structure of the matrix and we wonder if this structure can be implemented</div><div>in Petsc.</div><div><br></div><div>By special structure we mean the following. An entry in the linear system</div><div>is of the form (i, j) with 1<=i<=1000 and 1<=j<=N   with N = 3.10^5.</div><div>The node (i , j) is adjacent to all the nodes (i' , j) and thus they make a block</div><div>diagonal entry. But the node (i , j) is also adjacent to some nodes (i , j')</div><div>[About 6 such nodes, but it varies].</div></div></blockquote><div><br></div><div>I do not understand this explanation "(i, j) with 1<=i<=1000 and 1<=j<=N".</div><div>Your linear system is rectangular of size (1000, N)? Do you mean instead that each</div><div>entry in the linear system (i, j) is a 1000x1000 block?</div><div><br></div><div>We might have something that can help. If the structure of each entry is the same, we</div><div>have this class</div><div><br></div><div>  <a href="https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MATKAIJ.html">https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MATKAIJ.html</a></div><div><br></div><div>This would be writing the system as a Kronecker product</div><div><br></div><div>  A \otimes T</div><div><br></div><div>where T is your 1000x1000 matrix. We have only run this for moderate size T, say 16x16,</div><div>so further optimizations might be necessary, but it is a place to start.</div><div><br></div><div>Is it possible to write your system in this way?</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 dir="ltr"><div>Would there be a way to exploit this special structure in Petsc? I think</div><div>this should be fairly common and significant speedup could be obtained.</div><div><br></div><div>Best,</div><div><br></div><div>  Mathieu</div></div>
</blockquote></div><br clear="all"><div><br></div>-- <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>