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 <[email protected]> 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 <[email protected]> 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 <[email protected]> wrote:
>
> > On 27 September 2013 21:23, Barry Smith <[email protected]> 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
