Barry,
Thanks, that is essentially what I have done. I just modified it so
that before the call to KSPSolve it checks to make sure that max_it
for mglevels->smoothd is greater than zero. However, I still can't
think of a scenario where you might want to call KSPSolve with max_it
== 0 where you would want it to setup the PC, zero the solution
vector, etc., so I still think having KSPSolve exit immediately if
max_it==0 might be better.
John
On Mon, Aug 13, 2012 at 6:53 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
>
> 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
>>>>>
>>>>
>>>
>
