[petsc-dev] Branches with additional features

[Pierre Jolivet]

Also, is there something like MatDensePlaceArray in a branch somewhere, 
or is it possible to wrap a MatDense around a Vec?
I'm calling MatMatMult from on external Krylov solver (which does its 
own memory allocation), so at each iteration I either need to create new 
Mats with MatMPIDenseSetPreallocation, or use buffer Mats and use 
MatDenseGetArray/MatDenseRestoreArray. Both approaches have drawbacks I 
believe (the former needs MatAssemblyBegin/MatAssemblyEnd which does a 
global reduction, the latter is a waste of memory and needs copying 
from/to the memory allocated by PETSc) that would be avoided with 
something like MatDensePlaceArray.

Thanks.

On Mon, 22 May 2017 08:11:30 +0200, Pierre Jolivet wrote:
> On Sun, 21 May 2017 17:10:30 -0500, Barry Smith wrote:
>>> On May 21, 2017, at 3:17 PM, Pierre Jolivet 
>>> <pierre.jolivet at enseeiht.fr> wrote:
>>>
>>> Hello,
>>> I’m wondering if any branch has one or more of the following 
>>> features (or if this is in the pipeline somehow):
>>> 1) MKL PARDISO support for MPIAIJ
>>> 2) MKL PARDISO support for MPIBAIJ using iparm(37) 
>>> https://software.intel.com/en-us/mkl-developer-reference-fortran-pardiso-iparm-parameter#IPARM37
>>
>>   From:
>>
>>
>> 
>> https://software.intel.com/en-us/mkl-developer-reference-fortran-sparse-blas-bsr-matrix-storage-format#GUID-9BD4DADE-6059-4042-BC2A-F41048A47953
>>
>> values
>> A real array that contains the elements of the non-zero blocks of a
>> sparse matrix. The elements are stored block by block in row-major
>> order. A non-zero block is the block that contains at least one
>> non-zero element. All elements of non-zero blocks are stored, even 
>> if
>> some of them is equal to zero. Within each non-zero block the 
>> elements
>> are stored in column major order in the case of the one-based
>> indexing, and in row major order in the case of the zero-based
>> indexing.
>>
>> columns
>> Element i of the integer array columns is the number of the column 
>> in
>> the block matrix that contains the i-th non-zero block.
>>
>> rowIndex
>> Element j of this integer array gives the index of the element in 
>> the
>> columns array that is first non-zero block in a row j of the block
>> matrix.
>>
>> Unfortunately PETSc always stores "Within each non-zero block the
>> elements are stored in column major order" even though this it uses
>> zero based indexing. This means the matrix storage would need to be
>> monkeyed with before using MKL PARDISO. In order to reduce memory
>> usage and excess copying I would make copies of the columns and
>> rowIndex integers and shift them to one based indices and then you 
>> can
>> pass these to MKL Pardiso using the BCSR format; this is
>> straightforward for SeqBAIJ matrices. Like Matt said I don't know if
>> MKL PARDISO is appropriate for MPI matrices.
>
> OK, I see. That is true that MKL storage layout is somehow confusing
> as it depends on the numbering used.
> Concerning PARDISO, I got confused because Google only redirects to
> 
> http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MATSOLVERMKL_PARDISO.html#MATSOLVERMKL_PARDISO
> There is no such page as
> 
> http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MATSOLVERMKL_CPARDISO.html#MATSOLVERMKL_CPARDISO
> But MKL_CPARDISO is indeed defined there
> 
> http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatSolverPackage.html#MatSolverPackage,
> thanks for the pointer Stefano.
>
>>
>>> 3) MatMatMult_MPIBAIJ_MPIDense
>>
>>    We don't have any of the MatMatMult_XXXBAIJ_MPIDense() multiple
>> routines but they should be pretty easy to add, just copy the
>> corresponding AIJ routines and make the appropriate adjustments to
>> take into account the blocks.
>
> Yes, I've started copying most of the MatMatMult_XXXAIJ_XXXDense()
> routines. I was just wondering if someone already did that or tried 
> to
> use specialized libraries like libxsmm, but it looks like I'll have 
> to
> do it myself.
>
> Thanks for your inputs!
>
>>>
>>> Thanks in advance,
>>> Pierre