On Thu, Jul 5, 2018 at 8:32 AM, Mark Adams <[email protected]> wrote: > > > On Wed, Jul 4, 2018 at 1:09 PM Smith, Barry F. <[email protected]> wrote: > >> >> >> > On Jul 4, 2018, at 9:40 AM, Mark Adams <[email protected]> wrote: >> > >> > "-mat_seqaij_type seqaijmkl" just worked. >> >> Mark, >> >> Please clarify. Does this mean you can use -mat_seqaij_type >> seqaijmkl to satisfy all your needs now without changing any code? > > > Apparently yes. I timed the code and it took much longer with > NUM_OMP_THREADS > 1 (not sure what that is about but not my problem now), > so that tells me that threads were activated. > > I want to package up this simple approach and my branch that works with > -mat_type > aijmkl , and let the user analyze each. (So I am not worrying about > performance, I just want correctness). >
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. 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. 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. --Richard > >> >> >> Barry >> >> >> > On Wed, Jul 4, 2018 at 9:44 AM Mark Adams <[email protected]> wrote: >> > >> > >> > On Wed, Jul 4, 2018 at 3:01 AM Richard Tran Mills <[email protected]> >> 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 <[email protected]> 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 >> > >> >>
