On Sep 8, 2013, at 1:44 PM, Barry Smith <[email protected]> wrote: > > We should use this as an opportunity to fix a flaw in the original PC/MG and > Mat API. Note that PCMG uses Mat and Vec methods for everything EXCEPT this > one case > > ierr = PCMGSetResidual(pc,i,PCMGResidualDefault,mat);CHKERRQ(ierr); > > where the user can provide a matrix kernel as a function pointer to MG and MG > has a default version. Of course, no one actually changes from the default. > > We should add MatResidual() to the basic Mat methods (note this is the same > as MatMultSub() if we had such a thing) then all the many places in PETSc > where we have MatMult(); MatAXPY/AYPX(); can be replace with that one call > and PCMGSetResidual() and PCMGResidualDefault() would just disappear. We > could have MatResidual_Default() that did the usual MatMult followed by > VecAYPX() but each Mat implementation can implement an optimized one if they > want. >
So I should add a residual method to Mat, get rid of PCMGResidualDefault, make a MatResidualDefault, and have SetUp_MatAIJ set this default. I guess PCMGSetResidual would just set the residual for Mat to keep the same interface or should we just get rid of it and make a MatSetResidual? This used: 20:51 ~/Codes$ global -r PCMGSetResidual petsc/include/petscpc.h petsc/src/contrib/keyes/ex15.c petsc/src/ksp/pc/examples/tests/ex5.c petsc/src/ksp/pc/impls/mg/ftn-custom/zmgfuncf.c petsc/src/ksp/pc/impls/mg/mg.c petsc/src/ksp/pc/impls/ml/ml.c petsc/src/snes/examples/tests/ex11.c So I dive into these places and fix this up. And then we get rid of mglevels[l]->residual and replace it with mglevels[l]->A->residual. Is this what we have in mind? Mark > ---- > > After this change the SOR residual "optimization" becomes completely the > internal business of any particular matrix implementation. There would > naturally be both a MatResidual_SeqAIJ(), and MatResidual_SeqAIJ_inode(). Now > if MatSOR_SeqAIJ (and the inode version) would stash the partial row sums, > the vector it was applied to and the "state" of the matrix and the vector > then MatResidual_SeqAIJ() would check if the matrix and vector states match > the previous values and do the 1/2 version if they match, otherwise do the > standard residual computation. This way supports the optimization > automatically without the user or the MG code needing to have any special > checks or functions to call to use the optimization. > > Barry > > > On Sep 8, 2013, at 11:16 AM, "Mark F. Adams" <[email protected]> wrote: > >>>> >>>> 3) How are we going to intercept the residual call? We only use this for >>>> MG so we could modify the V-cycle to call a special residual method or >>>> check the matrix to see if has the upper part stashed… >>> >>> MGSetComputeResidual() already exists so this is how you would prorovide >>> the "special residual method" >>> >> >> Currently we have in PCSetUp_MG: >> >> ierr = PCMGSetResidual(pc,i,PCMGResidualDefault,mat);CHKERRQ(ierr); >> >> with >> >> PetscErrorCode PCMGSetResidual(PC pc,PetscInt l,PetscErrorCode >> (*residual)(Mat,Vec,Vec,Vec),Mat mat) >> { >> …. >> if (residual) mglevels[l]->residual = residual; >> if (!mglevels[l]->residual) mglevels[l]->residual = PCMGResidualDefault; >> …. >> >> PetscErrorCode PCMGResidualDefault(Mat mat,Vec b,Vec x,Vec r) >> { >> ... >> ierr = MatMult(mat,x,r);CHKERRQ(ierr); >> ierr = VecAYPX(r,-1.0,b);CHKERRQ(ierr); >> …. >> >> I propose: >> >> 1) add a flag, 'MatResidualType res_type;', to Mat with: >> >> typedef enum {MAT_RES_NONE, MAT_RES_SOR} MatResidualType; >> >> 2) Have MatSOR_SeqAIJ, etc., set this to MAT_RES_SOR if it cached data. >> >> 3) Have PCMGResidualDefault switch on mat->res_type, to call old code for >> MAT_RES_NONE and the new SOR code for MAT_RES_SOR. THe new code could be >> MatSORResidual(mat,b,x,r) and would have to switch between the blocked >> (inode.c) and regular (aij.c) version. >> >> This put the responsibility of making sure there is a valid cache on the >> user, in that the user must use this residual with the results of the >> smoother call. >> >> Mark >
