<div dir="ltr"><div><div>Hi Barry,<br></div>thank you for your feedback. I could be wrong but it seems to me that <span style="background:none repeat scroll 0% 0% yellow" class="">MatSetValuesBlocked</span> has to be called passing a 1D array, as you said.<br>
<br></div>In fact, From <a href="http://www">http://www</a>.<span style="background:none repeat scroll 0% 0% yellow" class="">mcs</span>.<span style="background:none repeat scroll 0% 0% yellow" class="">anl</span>.gov/<span style="background:none repeat scroll 0% 0% yellow" class="">petsc</span>/<span style="background:none repeat scroll 0% 0% yellow" class="">petsc</span>-current/docs/<span style="background:none repeat scroll 0% 0% yellow" class="">manualpages</span>/Mat/<span style="background:none repeat scroll 0% 0% yellow" class="">MatSetValuesBlocked</span>.html<br>
<div><pre> Suppose m=n=2 and block size(<span style="background:none repeat scroll 0% 0% yellow" class="">bs</span>) = 2 The array is
</pre>
<pre> 1 2 | 3 4
</pre>
<pre> 5 6 | 7 8
</pre>
<pre> - - - | - - -
</pre>
<pre> 9 10 | 11 12
</pre>
<pre> 13 14 | 15 16
</pre>
<pre> v[] should be passed in like
</pre>
<pre> v[] = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]
</pre>
<pre> If you are not using row oriented storage of v (that is you called <a href="http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatSetOption.html#MatSetOption"><span style="background:none repeat scroll 0% 0% yellow" class="">MatSetOption</span></a>(mat,MAT_ROW_ORIENTED,<a href="http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PETSC_FALSE.html#PETSC_FALSE"><span style="background:none repeat scroll 0% 0% yellow" class="">PETSC</span>_FALSE</a>)) then
</pre>
<pre> v[] = [1,5,9,13,2,6,10,14,3,7,11,15,4,8,12,16]
<br><br></pre>So, from above I can see that v[] is a 1D array. My problem or question is that my small block are stored in memory as<br></div><div> <br>
</div><div>A = [a11, a12 ; a21, a22] </div><div>B = [b11, b12 ; b21, b22] <br></div><div>etc.<br> (My actual blocks are 5x5 as indicated in my first e-mail but here I made them smaller to type less)<br><br></div>
<div><br></div><div>They are actually 2D arrays, i.e with two indices. Therefore before to pass them to <span style="background:none repeat scroll 0% 0% yellow" class="">MatSetValuesBlocked</span> I have to concatenate them in row oriented storage or column oriented storage, i.e<br>
<br></div><div>v = [a11, a12, b11, b12, a21, a22, b21, b22] <br><br></div><div>or <br><br></div><div>v = [a11, a21, a12, a22, b11, b21, b12, b22] <br></div><div><br>For each element in my mesh I construct N [5x5] blocks (N 2D arrays) which, once inserted in the PC matrix, will form the element contribution to the PC matrix. <font size="1"><font><b><span style="color:rgb(0,0,255)"></span></b></font></font>The 2D arrays are dense matrix and they are the output of the
multiplication between a dense matrix and a matrix in <span style="background:none repeat scroll 0% 0% yellow" class="">CSR</span> which is done
internally to my <span style="background:none repeat scroll 0% 0% yellow" class="">Fortran</span> code. Thus I will loop over all elements owned by a processor, and each time I will construct the N 2D arrays of 1 element and insert them in the PC matrix.<br>
</div><div><br></div><div>I would like to know if I can avoid to concatenate them in row oriented storage or column oriented storage and pass them directly as 2D arrays. <br><br></div><div><br></div><div>Thus, if I am not wrong, it would be nice if the user could have the option to pass to <span style="background:none repeat scroll 0% 0% yellow" class="">MatSetValuesBlocked</span> each block as:<br>
<br></div><div>1- values[<span style="background:none repeat scroll 0% 0% yellow" class="">bs</span>][<span style="background:none repeat scroll 0% 0% yellow" class="">bs</span>][1] (e.g. A = [a11, a12 ; a21, a22] ) <br>
</div><div> values[<span style="background:none repeat scroll 0% 0% yellow" class="">bs</span>][<span style="background:none repeat scroll 0% 0% yellow" class="">bs</span>][2] (e.g. B = [b11, b12 ; b21, b22] )<br></div>
<div><br>plus the two vectors that complete the <span style="background:none repeat scroll 0% 0% yellow" class="">CSR</span> format, i.e. global block row and column indices (i and j indices). Maybe the last indices that define the blocks (i.e. [1] and [2] in the examples can be avoided)<br>
</div><div><br>or <br><br>2- values[<span style="background:none repeat scroll 0% 0% yellow" class="">nblockrows</span>][<span style="background:none repeat scroll 0% 0% yellow" class="">nblockcols</span>][<span style="background:none repeat scroll 0% 0% yellow" class="">bs</span>][<span style="background:none repeat scroll 0% 0% yellow" class="">bs</span>] <br>
</div><div><br></div><div>Thanks,<br></div><div><br></div></div><div class="gmail_extra"><br><br><div class="gmail_quote">On Fri, Sep 27, 2013 at 5:34 PM, Barry Smith <span dir="ltr"><<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><br>
I don't understand the (mild) rant below, the expected format of the values is suppose to be EXACTLY what one puts in the double *array for block CSR format. That is the values are (formally) a 1 dimensional array) containing each block of values (in column oriented storage, that is the first column of the values in the first block, followed by the next column in the first block, etc etc) in the same order as the j indices. Note that this means that in Fortran the values array could be thought of as double A(bs,bs,nnz) Note that this is exactly how Matteo has his stuff organized so "it should just work". If it does not then please submit a bug report with a very small matrix that reproduces the problem.<br>
<br>
This function is exactly what it says it is, it takes in exactly a matrix in block CSR format.<br>
<div class="im"><br>
> wondering if a layout like values[nblockrows][nblockcols][bs][bs]<br>
> would not be much more convenient for users.<br>
<br>
</div> I don't understand this. It is a dense matrix as big as the entire sparse matrix.<br>
<span class="HOEnZb"><font color="#888888"><br>
Barry<br>
</font></span><div class="HOEnZb"><div class="h5"><br>
<br>
On Sep 27, 2013, at 4:02 PM, Lisandro Dalcin <<a href="mailto:dalcinl@gmail.com">dalcinl@gmail.com</a>> wrote:<br>
<br>
> On 27 September 2013 21:23, Barry Smith <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>> wrote:<br>
>> Just use MatCreateMPIBAIJWithArrays() saves thinking and work.<br>
><br>
> BTW, Have you ever thought about how hard to use is this API if you<br>
> want to provide the matrix values? How are users supposed to fill the<br>
> "values" array for a matrix with bs>1 ? This API is certainly not<br>
> Fortran-friendly (and I think it not even C-friendly). Mateo had the<br>
> CSR indices as well as the block values, but I recommended to use the<br>
> CSR indices to perform proper matrix preallocation, then loop over<br>
> cells and call MatSetValuesBlocked(). Otherwise, he would need to<br>
> create an new intermediate array and fill it the right way for<br>
> MatCreateMPIBAIJWithArrays() to work as expected.<br>
><br>
> The root of the problem is how MatSetValuesBlocked() expects the array<br>
> of values to be layout in memory. While I understand that the current<br>
> convention make it compatible with MatSetValues(), I'm always<br>
> wondering if a layout like values[nblockrows][nblockcols][bs][bs]<br>
> would not be much more convenient for users.<br>
><br>
><br>
><br>
> --<br>
> Lisandro Dalcin<br>
> ---------------<br>
> CIMEC (UNL/CONICET)<br>
> Predio CONICET-Santa Fe<br>
> Colectora RN 168 Km 472, Paraje El Pozo<br>
> 3000 Santa Fe, Argentina<br>
> Tel: <a href="tel:%2B54-342-4511594" value="+543424511594">+54-342-4511594</a> (ext 1016)<br>
> Tel/Fax: <a href="tel:%2B54-342-4511169" value="+543424511169">+54-342-4511169</a><br>
<br>
</div></div></blockquote></div><br><br clear="all"><br>-- <br><div dir="ltr">Matteo<br></div>
</div>