[petsc-dev] sor smoothers

[Barry Smith]

On Sep 8, 2013, at 7:53 PM, "Mark F. Adams" <mfadams at lbl.gov> wrote:

> 
> On Sep 8, 2013, at 1:44 PM, Barry Smith <bsmith at mcs.anl.gov> 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" <mfadams at lbl.gov> 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
>> 
>