On Fri, Jul 6, 2018 at 2:26 PM, Mark Adams <mfad...@lbl.gov> wrote: > >> Glad to hear that you got something (at least sort of) working, Mark. Let >> me know if you want me to work with your users on understanding the >> performance they are seeing. When I was working for Intel, I probably spent >> more time on dealing with stupid thread affinity settings than anything >> else when "optimizing" code performance. >> >> > That is a good idea. I will loop you into the thread. > > >> I should point out that, even with all of the thread settings correct, >> the current PETSc code base might not get much out of OpenMP threading with >> MKL. I'm assuming since you are using GAMG that a significant portion of >> time is spent in MatPtAP() operations? I'm not sure how much use of the >> sequential matrix-matrix multiply code the parallel MatPtAP() code can >> achieve. >> > > The RAP or AP is in MKL right? I understand they have these sparse > matrix-matrix kernels. The get some reuse out of the RAP I think. I will > ask when I bring you in the thread (but I will wait for you to clarify what > is supported). >
My SeqAIJMKL routines currently implement MatMatMult(), MatMatMultNumeric(), MatTransposeMatMult(), MatPtAP(), and MatPtAPNumeric(). Notice that there are no MatMatMultSymbolic and MatPtAPSymbolic() routines in there. The problem is that, even though MKL internally computes a symbolic multiplication and then a numeric one, until version 18 update 2 (the latest MKL), there was only one interface provided: one that did the full (symbolic + numeric) multiplication, that is, that maps to MatMatMult(). I (and I think a few others) complained about this, so in version 18 update 2, they added a "two stage" interface. Unfortunately, it's not the two stages that I want: The first stage only determines the amount of memory needed but doesn't compute the sparsity pattern. I did add some "Numeric" routines that use the sparsity pattern from a previous full multiplication (so for the MAT_REUSE_MATRIX case), but I've recently learned that what I've done isn't actually safe because that causes a memory leak in MKL. Frustrating. Support for a MatPtAP() operation was added very recently in MKL, but it is only for a symmetric matrix A. I check for MatIsSymmetricKnown(), and call the MKL version only if this check returns true; otherwise, I just call MatPtAP_SeqAIJ_SeqAIJ(). --Richard > >> There are some issues with mapping between the MKL interfaces for this >> and the PETSc ones, e.g., there is currently no equivalent of >> MatMultSymbolic() in MKL. I've had some discussions with the MKL team about >> this, but the current model they have isn't one that I can map all of the >> PETSc matrix-matrix interface onto. Hopefully I can get them to change this. >> > > OK, I don't think it will be a problem in the near term ... and I don't > understand why there is no equivalent of MatMultSymbolic (setting up > communication maps?). How does if function without this. Please clarify! > > Thanks, > > >> >> --Richard >> >> >>> >>>> >>>> >>>> Barry >>>> >>>> >>>> > On Wed, Jul 4, 2018 at 9:44 AM Mark Adams <mfad...@lbl.gov> wrote: >>>> > >>>> > >>>> > On Wed, Jul 4, 2018 at 3:01 AM Richard Tran Mills <rtmi...@anl.gov> >>>> wrote: >>>> > Hi Mark, >>>> > >>>> > I'd never looked at the code for MatCreateMPIAIJWithSeqAIJ(), but it >>>> looks like if you are using MPIAIJ matrices but you've set >>>> "-mat_seqaij_type seqaijmkl", >>>> > >>>> > No, I am using -mat_type aijmkl >>>> > >>>> > then the "diagonal" portion of the MPIAIJ matrix ends up as a >>>> SEQAIJMKL instance, but the off-diagonal portion ends up as a plain SEQAIJ >>>> when it is created with MatCreateSeqAIJWithArrays(). (I have not stepped >>>> through this with a debugger to verify that this is actually what happens.) >>>> So, if we want the off-diagonal matrix to also be a SEQAIJMKL, then we >>>> either need to add some logic after the MatCreateSeqAIJWithArrays() call to >>>> check to see if a different subtype of SEQAIJ is being used and then >>>> convert it to the correct type, or we need a variant of >>>> MatCreateSeqAIJWithArrays() that will honor the subtype of SEQAIJ that we >>>> are using. Barry, since you added the concept of subtyping for SEQAIJ, what >>>> do you prefer? >>>> > >>>> > No preference. >>>> > >>>> > I note that in some (most?) cases, I think we are OK or even better >>>> off with the off-diagonal matrix in the MPIAIJ not being a SEQAIJMKL one, >>>> because the PETSc routines can used the "compressed row" stuff to speed up >>>> handling the many zero rows that often end up in the off-diagonal, whereas >>>> the MKL routines cannot. >>>> > >>>> > >>>> > OK , I'll keep that in mind. I pushed a branch with fixes to get my >>>> gamg tests working (ksp ex56 primarily). I'm going to make a new branch and >>>> try "-mat_seqaij_type seqaijmkl". >>>> > >>>> > Thanks, >>>> > Mark >>>> > >>>> > --Richard >>>> > >>>> > On Tue, Jul 3, 2018 at 4:31 PM, Mark Adams <mfad...@lbl.gov> wrote: >>>> > I have having to fix AIJ methods that don't get the type from the >>>> args to set created matrix type, so as to keep the MKL typing. >>>> > >>>> > I am not sure how to fix this because I don't know to say 'make MPI >>>> from SEQ' in this method >>>> > >>>> > PetscErrorCode MatCreateMPIAIJWithSeqAIJ(MPI_Comm comm,Mat A,Mat >>>> B,const PetscInt garray[],Mat *mat) >>>> > { >>>> > .... >>>> > ierr = MatCreate(comm,mat);CHKERRQ(ierr); >>>> > ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); >>>> > .... >>>> > ierr = MatSetSizes(*mat,m,n,PETSC_DECIDE,N);CHKERRQ(ierr); >>>> > ierr = MatSetType(*mat,MATMPIAIJ);CHKERRQ(ierr); >>>> > >>>> > How should I do this? Other places in the code I get the type from >>>> input (MPI) args and set created matrix types accordingly. I could just put >>>> a switch with an error for any new types that might come along in the >>>> future ... >>>> > >>>> > Mark >>>> > >>>> >>>> >>