<div dir="ltr"><div dir="ltr">On Fri, Mar 17, 2023 at 2:23 PM Berger Clement <<a href="mailto:clement.berger@ens-lyon.fr">clement.berger@ens-lyon.fr</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 style="font-size:10pt;font-family:Verdana,Geneva,sans-serif">
<p>My issue is that it seems to improperly with some step of my process, the solve step doesn't provide the same result depending on the number of processors I use. I manually tried to multiply one the matrices I defined as a nest against a vector, and the result is not the same with e.g. 1 and 3 processors. That's why I tried the toy program I wrote in the first place, which highlights the misplacement of elements.</p></div></blockquote><div>Ah, now I think I understand.</div><div><br></div><div>  The PETSC_DECIDE arguments for sizes change with a different number of processes. You can put in har numbers if you want.</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 style="font-size:10pt;font-family:Verdana,Geneva,sans-serif">
<div>---<br>
<div style="margin:0px;padding:0px;font-family:monospace">Clément BERGER<br> ENS de Lyon</div>
</div>
<p><br></p>
<p>Le 2023-03-17 19:14, Barry Smith a écrit :</p>
<blockquote type="cite" style="padding:0px 0.4em;border-left:2px solid rgb(16,16,255);margin:0px"><br>
<div>   This sounds  like a fine use of MATNEST. Now back to the original question</div>
<div> </div>
<div>
<blockquote type="cite" style="padding:0px 0.4em;border-left:2px solid rgb(16,16,255);margin:0px">
<div style="font-size:10pt;font-family:Verdana,Geneva,sans-serif">
<blockquote style="padding:0px 0.4em;border-left-color:rgb(16,16,255);margin:0px">
<blockquote style="padding:0px 0.4em;border-left-color:rgb(16,16,255);margin:0px">
<div style="font-size:10pt">
<blockquote style="padding:0px 0.4em;border-left:2px solid rgb(16,16,255);margin:0px">
<blockquote style="padding:0px 0.4em;border-left:2px solid rgb(16,16,255);margin:0px">
<div style="font-size:10pt">
<p>I want to construct a matrix by blocs, each block having different sizes and partially stored by multiple processors. If I am not mistaken, the right way to do so is by using the MATNEST type. However, the following code</p>
<p>Call MatCreateConstantDiagonal(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,4,4,2.0E0_wp,A,ierr)<br>Call MatCreateConstantDiagonal(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,4,4,1.0E0_wp,B,ierr)<br>Call MatCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL_INTEGER,2,PETSC_NULL_INTEGER,(/A,PETSC_NULL_MAT,PETSC_NULL_MAT,B/),C,ierr)</p>
<p>does not generate the same matrix depending on the number of processors. It seems that it starts by everything owned by the first proc for A and B, then goes on to the second proc and so on (I hope I am being clear).</p>
<p>Is it possible to change that ?</p>
</div>
</blockquote>
</blockquote>
</div>
</blockquote>
</blockquote>
</div>
</blockquote>
  If I understand correctly it is behaving as expected. It is the same matrix on 1 and 2 MPI processes, the only difference is the ordering of the rows and columns. </div>
<div> </div>
<div>  Both matrix blocks are split among the two MPI processes. This is how MATNEST works and likely what you want in practice. </div>
<div><br>
<blockquote type="cite" style="padding:0px 0.4em;border-left:2px solid rgb(16,16,255);margin:0px">
<div>On Mar 17, 2023, at 1:19 PM, Berger Clement <<a href="mailto:clement.berger@ens-lyon.fr" target="_blank">clement.berger@ens-lyon.fr</a>> wrote:</div>
<br>
<div>
<div style="font-size:10pt;font-family:Verdana,Geneva,sans-serif">
<p>I have a matrix with four different blocks (2rows - 2columns). The block sizes differ from one another, because they correspond to a different physical variable. One of the block has the particularity that it has to be updated at each iteration. This update is performed by replacing it with a product of multiple matrices that depend on the result of the previous iteration. Note that these intermediate matrices are not square (because they also correspond to other types of variables), and that they must be completely refilled by hand (i.e. they are not the result of some simple linear operations). Finally, I use this final block matrix to solve multiple linear systems (with different righthand sides), so for now I use MUMPS as only the first solve takes time (but I might change it).</p>
<p>Considering this setting, I created each type of variable separately, filled the different matrices, and created different nests of vectors / matrices for my operations. When the time comes to use KSPSolve, I use MatConvert on my matrix to get a MATAIJ compatible with MUMPS, I also copy the few vector data I need from my nests in a regular Vector, I solve, I get back my data in my nest and carry on with the operations needed for my updates.</p>
<p>Is that clear ? I don't know if I provided too many or not enough details.</p>
<p>Thank you</p>
<div>---<br>
<div style="margin:0px;padding:0px;font-family:monospace">Clément BERGER<br> ENS de Lyon</div>
</div>
<p><br></p>
<p>Le 2023-03-17 17:34, Barry Smith a écrit :</p>
<blockquote style="padding:0px 0.4em;border-left:2px solid rgb(16,16,255);margin:0px">
<div> </div>
   Perhaps if you provide a brief summary of what you would like to do and we may have ideas on how to achieve it. 
<div> </div>
<div>   Barry</div>
<div> </div>
<div>Note: that MATNEST does require that all matrices live on all the MPI processes within the original communicator. That is if the original communicator has ranks 0,1, and 2 you cannot have a matrix inside MATNEST that only lives on ranks 1,2 but you could have it have 0 rows on rank zero so effectively it lives only on rank 1 and 2 (though its communicator is all three ranks).<br>
<div><br>
<blockquote style="padding:0px 0.4em;border-left:2px solid rgb(16,16,255);margin:0px">
<div>On Mar 17, 2023, at 12:14 PM, Berger Clement <<a href="mailto:clement.berger@ens-lyon.fr" target="_blank">clement.berger@ens-lyon.fr</a>> wrote:</div>
<br>
<div>
<div style="font-size:10pt;font-family:Verdana,Geneva,sans-serif">
<p>It would be possible in the case I showed you but in mine that would actually be quite complicated, isn't there any other workaround ? I precise that I am not entitled to utilizing the MATNEST format, it's just that I think the other ones wouldn't work.</p>
<div>---<br>
<div style="margin:0px;padding:0px;font-family:monospace">Clément BERGER<br> ENS de Lyon</div>
</div>
<p><br></p>
<p>Le 2023-03-17 15:48, Barry Smith a écrit :</p>
<blockquote style="padding:0px 0.4em;border-left:2px solid rgb(16,16,255);margin:0px">
<div> </div>
   You may be able to mimic what you want by not using PETSC_DECIDE but instead computing up front how many rows of each matrix you want stored on each MPI process. You can use 0 for on certain MPI processes for certain matrices if you don't want any rows of that particular matrix stored on that particular MPI process.
<div> </div>
<div>  Barry</div>
<div><br>
<div><br>
<blockquote style="padding:0px 0.4em;border-left:2px solid rgb(16,16,255);margin:0px">
<div>On Mar 17, 2023, at 10:10 AM, Berger Clement <<a href="mailto:clement.berger@ens-lyon.fr" target="_blank">clement.berger@ens-lyon.fr</a>> wrote:</div>
<div>
<div style="font-size:10pt;font-family:Verdana,Geneva,sans-serif">
<p>Dear all,</p>
<p>I want to construct a matrix by blocs, each block having different sizes and partially stored by multiple processors. If I am not mistaken, the right way to do so is by using the MATNEST type. However, the following code</p>
<p>Call MatCreateConstantDiagonal(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,4,4,2.0E0_wp,A,ierr)<br>Call MatCreateConstantDiagonal(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,4,4,1.0E0_wp,B,ierr)<br>Call MatCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL_INTEGER,2,PETSC_NULL_INTEGER,(/A,PETSC_NULL_MAT,PETSC_NULL_MAT,B/),C,ierr)</p>
<p>does not generate the same matrix depending on the number of processors. It seems that it starts by everything owned by the first proc for A and B, then goes on to the second proc and so on (I hope I am being clear).</p>
<p>Is it possible to change that ?</p>
<p>Note that I am coding in fortran if that has ay consequence.</p>
<p>Thank you,</p>
<p>Sincerely,</p>
<div>-- <br>
<div style="margin:0px;padding:0px;font-family:monospace">Clément BERGER<br> ENS de Lyon</div>
</div>
</div>
</div>
</blockquote>
</div>
</div>
</blockquote>
</div>
</div>
</blockquote>
</div>
</div>
</blockquote>
</div>
</div>
</blockquote>
</div>
</blockquote>
</div>
</blockquote></div><br clear="all"><div><br></div><span class="gmail_signature_prefix">-- </span><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>