[petsc-users] MatMPIBAIJSetPreallocationCSR i and j indices

Fri Sep 27 21:34:57 CDT 2013

  Mattero

    Multi-dimensional arrays in Fortran are are a sham. That is, a multi-dimensional array in fortran is simply a one dimensional array with a little additional information. 

   If you have a three dimensional array in Fortran that is 

    1 3    
    2 4 
          5 7
          6 8
                 9 11
                 10 12 

   where above I have indicated

    A(:,:,1)
              A(:,:,2)
                         A(:,:,3)

     it is in actual memory stored in consecutive memory locations as 

      1 2 3 4 5 6 7 8 9 10 11 12 

     this is called column major format  because the columns are listed in consecutive memory before rows.  

    When you call MatCreateMPIBAIJWithArrays() it does not expect the values in row-major order, it expects them in column major order per block, which, as I understand it, is exactly what you have already computed.

    Please just try what I suggest and see what happens. Send all the output if it goes wrong, for a tiny problem. The beautiful thing which I love about programming is no one punishes you for trying something, even it if it is wrong. As I said before it should just work, if it doesn't then let us know EXACTLY what goes wrong and we'll see if there is a bug on our end or a misunderstanding. We are talking too much and not trying enough.

   Barry

On Sep 27, 2013, at 8:01 PM, Matteo Parsani <parsani.matteo at gmail.com> wrote:

> 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(mat,MAT_ROW_ORIENTED,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