John,
You are right we do not have completely optimized forms of all variants of
the multigrid algorithms.
You could make a modification to the routine
#undef __FUNCT__
#define __FUNCT__ "PCMGMCycle_Private"
PetscErrorCode PCMGMCycle_Private(PC pc,PC_MG_Levels
**mglevelsin,PCRichardsonConvergedReason *reason)
{
PC_MG *mg = (PC_MG*)pc->data;
PC_MG_Levels *mgc,*mglevels = *mglevelsin;
PetscErrorCode ierr;
PetscInt cycles = (mglevels->level == 1) ? 1 : (PetscInt)
mglevels->cycles;
PetscFunctionBegin;
if (mglevels->eventsmoothsolve) {ierr =
PetscLogEventBegin(mglevels->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);}
ierr = KSPSolve(mglevels->smoothd,mglevels->b,mglevels->x);CHKERRQ(ierr); /*
pre-smooth */
if (mglevels->eventsmoothsolve) {ierr =
PetscLogEventEnd(mglevels->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);}
if (mglevels->level) { /* not the coarsest grid */
if (mglevels->eventresidual) {ierr =
PetscLogEventBegin(mglevels->eventresidual,0,0,0,0);CHKERRQ(ierr);}
ierr =
(*mglevels->residual)(mglevels->A,mglevels->b,mglevels->x,mglevels->r);CHKERRQ(ierr);
if (mglevels->eventresidual) {ierr =
PetscLogEventEnd(mglevels->eventresidual,0,0,0,0);CHKERRQ(ierr);}
/* if on finest level and have convergence criteria set */
if (mglevels->level == mglevels->levels-1 && mg->ttol && reason) {
PetscReal rnorm;
ierr = VecNorm(mglevels->r,NORM_2,&rnorm);CHKERRQ(ierr);
if (rnorm <= mg->ttol) {
if (rnorm < mg->abstol) {
*reason = PCRICHARDSON_CONVERGED_ATOL;
ierr = PetscInfo2(pc,"Linear solver has converged. Residual norm %G
is less than absolute tolerance %G\n",rnorm,mg->abstol);CHKERRQ(ierr);
} else {
*reason = PCRICHARDSON_CONVERGED_RTOL;
ierr = PetscInfo2(pc,"Linear solver has converged. Residual norm %G
is less than relative tolerance times initial residual norm
%G\n",rnorm,mg->ttol);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
}
mgc = *(mglevelsin - 1);
if (mglevels->eventinterprestrict) {ierr =
PetscLogEventBegin(mglevels->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);}
ierr = MatRestrict(mglevels->restrct,mglevels->r,mgc->b);CHKERRQ(ierr);
if (mglevels->eventinterprestrict) {ierr =
PetscLogEventEnd(mglevels->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);}
ierr = VecSet(mgc->x,0.0);CHKERRQ(ierr);
while (cycles--) {
ierr = PCMGMCycle_Private(pc,mglevelsin-1,reason);CHKERRQ(ierr);
}
if (mglevels->eventinterprestrict) {ierr =
PetscLogEventBegin(mglevels->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);}
ierr =
MatInterpolateAdd(mglevels->interpolate,mgc->x,mglevels->x,mglevels->x);CHKERRQ(ierr);
if (mglevels->eventinterprestrict) {ierr =
PetscLogEventEnd(mglevels->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);}
if (mglevels->eventsmoothsolve) {ierr =
PetscLogEventBegin(mglevels->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);}
ierr = KSPSolve(mglevels->smoothu,mglevels->b,mglevels->x);CHKERRQ(ierr);
/* post smooth */
if (mglevels->eventsmoothsolve) {ierr =
PetscLogEventEnd(mglevels->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);}
}
PetscFunctionReturn(0);
}
to remove the unneeded computations.
Barry
On Aug 13, 2012, at 3:39 PM, John Fettig <john.fettig at gmail.com> wrote:
> What if you wanted to do a full cycle with no pre-smooths instead of a
> v-cycle?
>
> John
>
> On Mon, Aug 13, 2012 at 4:34 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
>>
>> #undef __FUNCT__
>> #define __FUNCT__ "PCMGKCycle_Private"
>> PetscErrorCode PCMGKCycle_Private(PC pc,PC_MG_Levels **mglevels)
>> {
>> PetscErrorCode ierr;
>> PetscInt i,l = mglevels[0]->levels;
>>
>> PetscFunctionBegin;
>> /* restrict the RHS through all levels to coarsest. */
>> for (i=l-1; i>0; i--){
>> if (mglevels[i]->eventinterprestrict) {ierr =
>> PetscLogEventBegin(mglevels[i]->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);}
>> ierr =
>> MatRestrict(mglevels[i]->restrct,mglevels[i]->b,mglevels[i-1]->b);CHKERRQ(ierr);
>> if (mglevels[i]->eventinterprestrict) {ierr =
>> PetscLogEventEnd(mglevels[i]->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);}
>> }
>>
>> /* work our way up through the levels */
>> ierr = VecSet(mglevels[0]->x,0.0);CHKERRQ(ierr);
>> for (i=0; i<l-1; i++) {
>> if (mglevels[i]->eventsmoothsolve) {ierr =
>> PetscLogEventBegin(mglevels[i]->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);}
>> ierr =
>> KSPSolve(mglevels[i]->smoothd,mglevels[i]->b,mglevels[i]->x);CHKERRQ(ierr);
>> if (mglevels[i]->eventsmoothsolve) {ierr =
>> PetscLogEventEnd(mglevels[i]->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);}
>> if (mglevels[i+1]->eventinterprestrict) {ierr =
>> PetscLogEventBegin(mglevels[i+1]->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);}
>> ierr =
>> MatInterpolate(mglevels[i+1]->interpolate,mglevels[i]->x,mglevels[i+1]->x);CHKERRQ(ierr);
>> if (mglevels[i+1]->eventinterprestrict) {ierr =
>> PetscLogEventEnd(mglevels[i+1]->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);}
>> }
>> if (mglevels[l-1]->eventsmoothsolve) {ierr =
>> PetscLogEventBegin(mglevels[l-1]->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);}
>> ierr =
>> KSPSolve(mglevels[l-1]->smoothd,mglevels[l-1]->b,mglevels[l-1]->x);CHKERRQ(ierr);
>> if (mglevels[l-1]->eventsmoothsolve) {ierr =
>> PetscLogEventEnd(mglevels[l-1]->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);}
>>
>> PetscFunctionReturn(0);
>> }
>>
>>
>> On Aug 13, 2012, at 3:19 PM, Jed Brown <jedbrown at mcs.anl.gov> wrote:
>>
>>> Shorthand for this is -pc_mg_type kaskade.
>>>
>>> On Mon, Aug 13, 2012 at 1:01 PM, John Fettig <john.fettig at gmail.com>
>>> wrote:
>>> Barry,
>>>
>>> Thank you for answering my question. I have another one for you: it
>>> seems the special case of zero pre-smooths is somewhat non-trivial.
>>> The best I can do is set the pre-smoother to Richardson with PCNONE
>>> and zero as max_its. However, if you aren't careful in setting
>>> KSPSetInitialGuessNonzero this can have unexpected results since the
>>> generic KSPSolve will clobber your solution before it even tries a
>>> convergence criteria (thus ruling out KSPPREONLY). It also does a
>>> couple of unnecessary residual calculations. Would it be reasonable to
>>> put a zero-iteration special case in KSPSolve so that if you don't
>>> want any iterations it doesn't actually do anything (no setup, no
>>> preconditioner, no residual, no scaling, etc.)?
>>>
>>> Thanks,
>>> John
>>>
>>> On Thu, Aug 9, 2012 at 6:37 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
>>>>
>>>> John,
>>>>
>>>> On Aug 9, 2012, at 9:50 AM, John Fettig <john.fettig at gmail.com> wrote:
>>>>
>>>>> I am a little confused about what Richardson means. If you use
>>>>> multiplicative V-cycle multigrid with Richardson KSP (and no
>>>>> convergence monitor), it sets the applyrichardson operator to
>>>>> PCApplyRichardson_MG, which appears to just run V-cycles until
>>>>> convergence.
>>>>
>>>> Yes, this is correct.
>>>>
>>>>> As far as I can tell, it doesn't ever update according
>>>>> to the documented
>>>>>
>>>>> x^{n+1} = x^{n} + scale*B(b - A x^{n})
>>>>>
>>>> In exact arithmetic it is actually "implicitly" doing exactly this
>>>> update. It is difficult to see why this is true generally (because B is
>>>> rather complicated for multigrid) but if you consider only two levels with
>>>> a direct solver on the coarse grid and SSOR as the pre and post smooth you
>>>> can write out the formulas and map back and forth between the two forms.
>>>> The reason for the PCApplyRichardson_ forms is because they are a bit more
>>>> efficient than separating out the action of B and then doing the update
>>>> as above.
>>>>
>>>>
>>>>> If on the other hand you use full MG, it does update according to the
>>>>> above formula. This also happens if you set a convergence monitor.
>>>>>
>>>>> I can see how multiplicative V-cycle with Richardson is simply using
>>>>> multigrid as a solver. What I don't understand is how full MG with
>>>>> Richardson is using multigrid as a solver, because it is using the
>>>>> update formula above in between cycles.. Shouldn't there be a
>>>>> applyrichardson for full multigrid as well? If not, why?
>>>>
>>>> I think there could be a applyRichardson for full multigrid but it
>>>> would be kind of complicated and would not benefit much because the amount
>>>> of work in a full multigrid step is much higher so the savings would be a
>>>> much lower percentage than with V cycle.
>>>>
>>>> Barry
>>>>
>>>>>
>>>>> Thanks,
>>>>> John
>>>>
>>>
>>
