[petsc-users] MatMPIBAIJSetPreallocationCSR i and j indices

[Jed Brown]

Barry Smith <bsmith at mcs.anl.gov> writes:

>    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) 

No, this is not how it works.  The ordering is

  values[nblockrows][bs][nblockcols][bs]

or, in Fortran

  values(bs,nblockcols,bs,nblockrows)

where nblockrows and nblockcols are the number of rows and columns being
inserted in this all to MatSetValuesBlocked.

this layout means that the entries are packed in exactly the same way as
if the block indices were blown up to scalar indices and then plain
MatSetValues was called.

>    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.
>
>    This function is exactly what it says it is, it takes in exactly a matrix in block CSR format.
>
>> wondering if a layout like values[nblockrows][nblockcols][bs][bs]
>> would not be much more convenient for users.
>
>   I don't understand this. It is a dense matrix as big as the entire sparse matrix.

Lisandro means the number of block rows and cols being inserted in the
call, not the dimension of the matrix.

>
>    Barry
>
>   
> On Sep 27, 2013, at 4:02 PM, Lisandro Dalcin <dalcinl at gmail.com> wrote:
>
>> On 27 September 2013 21:23, Barry Smith <bsmith at mcs.anl.gov> wrote:
>>> Just use MatCreateMPIBAIJWithArrays() saves thinking and work.
>> 
>> BTW, Have you ever thought about how hard to use is this API if you
>> want to provide the matrix values? How are users supposed to fill the
>> "values" array for a matrix with bs>1 ? This API is certainly not
>> Fortran-friendly (and I think it not even C-friendly). Mateo had the
>> CSR indices as well as the block values, but I recommended to use the
>> CSR indices to perform proper matrix preallocation, then loop over
>> cells and call MatSetValuesBlocked(). Otherwise, he would need to
>> create an new intermediate array and fill it the right way for
>> MatCreateMPIBAIJWithArrays() to work as expected.
>> 
>> The root of the problem is how MatSetValuesBlocked() expects the array
>> of values to be layout in memory. While I understand that the current
>> convention make it compatible with MatSetValues(), I'm always
>> wondering if a layout like values[nblockrows][nblockcols][bs][bs]
>> would not be much more convenient for users.
>> 
>> 
>> 
>> -- 
>> Lisandro Dalcin
>> ---------------
>> CIMEC (UNL/CONICET)
>> Predio CONICET-Santa Fe
>> Colectora RN 168 Km 472, Paraje El Pozo
>> 3000 Santa Fe, Argentina
>> Tel: +54-342-4511594 (ext 1016)
>> Tel/Fax: +54-342-4511169
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 835 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20130928/0e8fab4b/attachment-0001.pgp>