[petsc-users] MatMPIBAIJSetPreallocationCSR i and j indices

[Matteo Parsani]

Hi Barry,
thank you for your feedback. I could be wrong but it seems to me that
MatSetValuesBlocked has to be called passing a 1D array, as you said.

In fact, From http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages
/Mat/MatSetValuesBlocked.html

  Suppose m=n=2 and block size(bs) = 2 The array is

  1  2  | 3  4

  5  6  | 7  8

  - - - | - - -

  9  10 | 11 12

  13 14 | 15 16

  v[] should be passed in like

  v[] = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]

 If you are not using row oriented storage of v (that is you called
MatSetOption <http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatSetOption.html#MatSetOption>(mat,MAT_ROW_ORIENTED,PETSC_FALSE
<http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PETSC_FALSE.html#PETSC_FALSE>))
then

  v[] = [1,5,9,13,2,6,10,14,3,7,11,15,4,8,12,16]

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


A = [a11, a12 ;  a21, a22]
B = [b11, b12 ; b21, b22]
etc.
  (My actual blocks are 5x5 as indicated in my first e-mail but here I made
them smaller to type less)


They are actually 2D arrays, i.e with two indices. Therefore before to pass
them to MatSetValuesBlocked I have to concatenate them in row oriented
storage or column oriented storage, i.e

v = [a11, a12, b11, b12, a21, a22, b21, b22]

or

v = [a11, a21, a12, a22, b11, b21, b12, b22]

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. **The 2D arrays are dense matrix and they are the output of the
multiplication between a dense matrix and a matrix in CSR which is done
internally to my Fortran 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.

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.


Thus, if I am not wrong, it would be nice if the user could have the option
to pass to MatSetValuesBlocked each block as:

1- values[bs][bs][1]  (e.g. A = [a11, a12 ;  a21, a22] )
    values[bs][bs][2]  (e.g. B = [b11, b12 ; b21, b22] )

plus the two vectors that complete the CSR 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)

or

2- values[nblockrows][nblockcols][bs][bs]

Thanks,



On Fri, Sep 27, 2013 at 5:34 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:

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


-- 
Matteo
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20130927/acf6410f/attachment.html>