[petsc-dev] Branches with additional features

[Barry Smith]

> On May 28, 2017, at 2:51 PM, Pierre Jolivet <Pierre.Jolivet at enseeiht.fr> wrote:
> 
> 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.
> 

   We don't currently have a MatDensePlaceArray() but could add it if needed. First, let me check, does your matrix change size or is it fixed size (MatDensePlaceArray() would probably only be possible if the size remains the same)?

   Barry

> 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