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