[petsc-dev] Branches with additional features

[Pierre Jolivet]

> On 29 May 2017, at 06:49, Jed Brown <jed at jedbrown.org> wrote:
> 
> You can't get a callback from this "external Krylov solver"?

You mean, a callback to let the external solver use the memory allocated by PETSc for the Mats? Unfortunately, I need the underlying memory to be contiguous between all Mats, and I don't think it is possible with the current PETSc operators for matrix creation.

> I'm not aware of the hack/feature being in a branch, but note that
> MatMatMult_MPIDense_MPIDense calls MatConvert_MPIDense_Elemental which
> copies values into the distribution for which the implementation can be
> efficient.  The overhead of your other copy is trivial in comparison
> (except perhaps when you have very few columns).

I'll look into that, thanks.

Barry, the size of the matrices are not changing, except sometimes their number of columns (they are tall and skinny matrices) when the Krylov iterations are restarted, but I could just re-create a matrix then (instead of at every iteration).

Thanks.

> Pierre Jolivet <Pierre.Jolivet at enseeiht.fr> writes:
> 
>> 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