On Sep 8, 2013, at 7:53 PM, "Mark F. Adams" <[email protected]> wrote:
> > 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. > MatResidual_Default MatCreate_SeqAIJ would set this default. In fact all MatCreate_xxx would need to set the 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? There would be no MatSetResidual() in the same way there is not MatSetMult and no PCMGSetResidual() since the Mat now "knows" how to compute a residual. > > 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. No, it would be MatResidual(mglevels[i]->A, ….) just like MatMult(mglevels[i]->A,….) etc > > 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 >> >
