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