Well, I thought of AIJMKL keeping track of the state for which 
mkl_sparse_optimize() was called. If the matrix state changes, the next call to 
MatMult()

    is MatMult the only operation that needs this check, or do many Mat methods 
need this check?

Looks like it is needed for MatMult(), MatMatMul() and triangular solves:
https://software.intel.com/en-us/mkl-developer-reference-fortran-inspector-executor-sparse-blas-analysis-routines


    For parallel matrices are all the operations that need to do this kind of 
update collective over the Mat?

According to the list above, all operations are collective.


will detect that the current state does not match the state it was optimized 
for and hence trigger another optimization. This isn't too different from what 
we do for GPU stuff.

    We do this for cusparse

static PetscErrorCode MatMult_SeqAIJCUSPARSE(Mat A,Vec xx,Vec yy)
{
....
   PetscFunctionBegin;
   /* The line below is necessary due to the operations that modify the matrix 
on the CPU (axpy, scale, etc) */
   ierr = MatSeqAIJCUSPARSECopyToGPU(A);CHKERRQ(ierr);

    But for CUSP (deprecated) we do it at MatAssembly time

PetscErrorCode MatAssemblyEnd_SeqAIJCUSP(Mat A,MatAssemblyType mode)
{
   PetscErrorCode ierr;

   PetscFunctionBegin;
   ierr = MatAssemblyEnd_SeqAIJ(A,mode);CHKERRQ(ierr);
   ierr = MatCUSPCopyToGPU(A);CHKERRQ(ierr);


   as we do it for ViennaCL

PetscErrorCode MatAssemblyEnd_SeqAIJViennaCL(Mat A,MatAssemblyType mode)
{
   PetscErrorCode ierr;

   PetscFunctionBegin;
   ierr = MatAssemblyEnd_SeqAIJ(A,mode);CHKERRQ(ierr);
   ierr = MatViennaCLCopyToGPU(A);CHKERRQ(ierr);

   For SeqAIJPERM we do that extra processing of the data in MatAssemblyEnd and 
don't need to do anything if only the numerical values change.

   So we seem to have some inconsistencies in how we handle this which will 
lead to more errors.

ViennaCL should adopt the CUSPARSE model then (and CUSP be removed).

Best regards,
Karli

Reply via email to