<div dir="ltr"><div><div class="gmail_extra"><span style="font-family:courier new,monospace"></span><div class="gmail_quote"><span style="font-family:courier new,monospace">On 8 November 2013 22:22, Jed Brown <span dir="ltr"><<a href="mailto:jedbrown@mcs.anl.gov" target="_blank">jedbrown@mcs.anl.gov</a>></span> wrote:</span><br>

<span style="font-family:courier new,monospace"></span><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div><div><span style="font-family:courier new,monospace">
> I would like to create matrix G and vector R such that<br>
><br>
> G*P = R<br>
><br>
> where<br>
><br>
> G = /Gx\<br>
>     \Gy/<br>
> and<br>
><br>
> R = /Rx\<br>
>     \Ry/<br>
><br>
> and vector R is created from the DMComposite of the DAs of U and V.<br></span>
<span style="font-family:courier new,monospace"><br></span>
</div></div><span style="font-family:courier new,monospace">DMCreateGlobalVector() gives you the vector R.  Then use MatCreate() and<br>
MatSetSizes() where the local and global row sizes match R and the<br>
column sizes match P.  For your staggered grid, you should be able to<br>
preallocate simply by using MatSeqAIJSetPreallocation() and<br>
MatMPIAIJSetPreallocation() with a constant row size of 2 for the<br>
"diagonal block" and 1 in the off-diagonal block.  (You can also<br>
preallocate exactly to be more accurate with memory, but may as well<br>
know you are assembling the correct thing first.)<br></span>
<span style="font-family:courier new,monospace"><br>
You can call DMCompositeGetGlobalISs() to get index sets holding the<br>
global indices for the U and V rows of your matrix, or just compute them<br>
>From the ownership range and the sizes.  The column indices will come<br>
>From the DMDA that holds P.<br></span>
</blockquote></div><span style="font-family:courier new,monospace"><br></span><br><div><div><span style="font-family:courier new,monospace">Hi Jed,<br><br></span></div><span style="font-family:courier new,monospace">Thank you very much! That worked.<br>

<br>I also have a couple of general questions regarding this:<br><br>1. In this case, </span><span style="font-family:courier new,monospace"><span style="font-family:courier new,monospace">I followed your advice and created a single matrix G of type MPIAIJ and assembled it. H</span>ow do I decide if I should rather use Block Matrices (MPIBAIJ) or MatGetLocalSubMatrix to set it up?<br>

<br></span></div><span style="font-family:courier new,monospace">2. Various functions exist to help us access the correct indices when required (e.g. </span><span style="font-family:courier new,monospace">DMDAGetAO, DMCompositeGetGlobalISs, DMDAGetCorners, VecGetOwnershipRange, etc.) And sometimes, there are multiple ways by which we can obtain the same indices. In such cases, are some of these functions preferable to others (either because they are faster or need to allocate less memory)?<br>

<br></span></div><div class="gmail_extra"><span style="font-family:courier new,monospace">Thank you,<br></span></div><div class="gmail_extra"><span style="font-family:courier new,monospace">Anush<br></span></div></div></div>