On Sun, Jul 30, 2017 at 5:42 AM, Karl Rupp <r...@iue.tuwien.ac.at> wrote:
> >> > >> 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. > I think this sort of "lazy" approach to performing the matrix optimization is ultimately what we should do for the mkl_sparse_optimize() calls: We should probably do the optimization step inside the MatAssemblyEnd sequence, since that seems like a reasonable thing to tack on, but afterwards we should rely on this check for changed values to trigger another optimization. (Once we are sure that we've got everything that updates the matrix actually marking that it has changed the matrix state.) I don't like doing the sparse optimization after every call to, say, MatDuplicate, since it's mostly likely that the user is going to modify the values of the duplicated matrix before performing any operations with the new matrix -- but we can't guarantee that, so in the current model I do the optimization step. --Richard > > > > 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 >
