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. ---- 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
