[petsc-dev] ML Reuse redux

[John Fettig]

I've emailed this list about this previously, but I'm revisiting it
and thought I'd bounce this off you again.  When you are solving a
series of equations, where the nonzero pattern stays the same but the
values change, shouldn't it make sense to offer an option for AMG
solvers which changes the matrices used on the levels but does not
recompute the restriction/interpolating matrices?  I'm thinking
something like (to replace the if(pc->setupcalled) branch in ml.c in
PCSetUp_ML()):

if (pc->setupcalled){

      Nlevels = pc_ml->Nlevels;
      fine_level = Nlevels - 1;
      level = fine_level - 1;

      gridctx = pc_ml->gridctx;
      gridctx[fine_level].A = A;
      for (mllevel=1; mllevel<Nlevels; mllevel++){
          MatPtAP(gridctx[level+1].A, gridctx[level].P,
MAT_REUSE_MATRIX, 1.0, &gridctx[level].A);
          if (level > 0){
              ierr =
PCMGSetResidual(pc,level,PCMGDefaultResidual,gridctx[level].A);CHKERRQ(ierr);
          }
          ierr =
KSPSetOperators(gridctx[level].ksp,gridctx[level].A,gridctx[level].A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
          level--;
      }
      ierr = PCMGSetResidual(pc,fine_level,PCMGDefaultResidual,gridctx[fine_level].A);CHKERRQ(ierr);
      ierr = KSPSetOperators(gridctx[fine_level].ksp,gridctx[fine_level].A,gridctx[fine_level].A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
      /* setupcalled is set to 0 so that MG is setup from scratch */
      pc->setupcalled = 0;
      ierr = PCSetUp_MG(pc);CHKERRQ(ierr);
      PetscFunctionReturn(0);
} else {
/* .... rest of the original PCSetUp_ML goes here ... */
}


Of course it shouldn't be Pt.A.P but rather R.A.P, but you get the
idea (aside: how do you compute R.A.P. with petsc? this doesn't seem
to be implemented).  I find just by doing this, the setup cost for ML
is halved for subsequent solves, and the preconditioner works almost
as well as rebuilding from scratch (what the original code does). Is
there some reason to not do this (or at least make an option that
allows one to do this)?  If I use SAME_PRECONDITIONER instead, the
preconditioner isn't nearly as effective.

Regards,
John