[petsc-dev] sor smoothers

[Mark F. Adams]

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