I forgot that MatSetValuesBlocked() takes row oriented input even within the blocks. So if we simply call MatSetOption(mat,MAT_ROW_ORIENTED,PETSC_FALSE) on the matrix after it is created then the storage should work as I originally planned.
Barry I found this discussion values(bs,nblockcols,bs,nblockrows) very confusing. MatMPIBAIJSetPreallocationCSR() doesn't take a square block of values so I didn't see what you are talking about. You must be talking about each internal call to MatSetValuesBlocked_MPIBAIJ(B,1,&row,ncols,icols,svals,INSERT_VALUES) where it puts in one block row of values with each call. On Sep 28, 2013, at 9:12 PM, Jed Brown <[email protected]> wrote: > Barry Smith <[email protected]> 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 <[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
