[petsc-users] ILU and Block Gauss Seidel smoothing

[Gianluca Meneghello]

Thanks Barry.

What I'm doing at the moment in order to pass the ordering to the PC is

PCFactorSetMatOrderingType(PC,"FESOrder");

where the FESOrder(...) function computes the ordering and has been
registered in the main function with

MatOrderingRegister("FESOrder",0,"FESOrder",FESOrder);

The same is done for al the different orderings I'm using inside the
PCSHELL suggested by Jed.

My understanding is that the FESOrder(...) function is evaluated every
time the PCApply(...) function is called.

If there was a way to store the IS containing the ordering — or any
other array — into the shell context (this I know how to do that) and
the pass the IS directly to the preconditioner instead of using
PCFactorSetMatOrderingType would be the best solution for me.

Thanks

Gianluca


On 12 July 2011 20:15, Barry Smith <bsmith at mcs.anl.gov> wrote:
>
> On Jul 12, 2011, at 3:47 AM, Gianluca Meneghello wrote:
>
>> Jed, thanks for your answer.
>>
>> This is exactly what I was looking for (it took me some time to test it...).
>>
>> Let me ask you a couple of additional questions:
>>
>> 1) in this configuration, is there an alternative way to pass the
>> ordering to the preconditioner? Becuse my ordering depends on the grid
>> an not on the solution, it is constant during all the computations and
>> I have no advantage on calling the function that creates multiple
>> times (the pc is built and destroyed several times).
>
>   Not sure exactly what you want here. You can put the ordering information in the shell context and just keep it there forever and don't need to to recompute.
>
>    PCShellSet/GetContext()?
>>
>> 2) what's the best way to make the shell pc selectable from the option
>> database, so that it can be replaced with more standard ones?
>>
>
>        /* here the user can provide -pc_type shell or ilu or whatever */
>        KSPSetFromOptions(ksp);
>        KSPGetPC(ksp,&pc);
>        /* the following lines will only do something if the user selected the PCType of shell otherwise they are ignored */
>
>        PCShellSetContext(pc,....)
>        PCShellSetSetUp(pc,....)
>        PCShellSetApply(pc,.....)
>
>
>   Barry
>
>
>> Thanks
>>
>> Gianluca
>>
>> On 1 July 2011 14:43, Jed Brown <jed at 59a2.org> wrote:
>>> On Fri, Jul 1, 2011 at 07:11, Gianluca Meneghello <gianmail at gmail.com>
>>> wrote:
>>>>
>>>> I've tried to implement Jed's solution as it was looking faster to
>>>> try. This is what I came up with for two different orders, FENOrder
>>>> and FWNOrder respectively
>>>>
>>>> 283       KSP ksp;
>>>> 284       PC pc;
>>>> 285
>>>> 286       ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);  CHKERRQ(ierr);
>>>> 287       ierr = KSPSetOptionsPrefix(ksp,"mg_");    CHKERRQ(ierr);
>>>> 288       ierr = KSPSetOperators(ksp,J,J,SAME_NONZERO_PATTERN);
>>>> CHKERRQ(ierr);
>>>> 289       ierr = KSPSetFromOptions(ksp);    CHKERRQ(ierr);
>>>> 290       ierr = KSPGetPC(ksp,&pc);         CHKERRQ(ierr);
>>>> 291       ierr = PCFactorSetFill(pc,1);     CHKERRQ(ierr);
>>>> 292       ierr = PCFactorSetMatOrderingType(pc,"FENOrder"); CHKERRQ(ierr);
>>>> 293       ierr = KSPSolve(ksp,b,x); CHKERRQ(ierr);
>>>> 294       ierr = KSPDestroy(&ksp);  CHKERRQ(ierr);
>>>> 295
>>>> 296       ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);  CHKERRQ(ierr);
>>>> 297       ierr = KSPSetOptionsPrefix(ksp,"mg_");    CHKERRQ(ierr);
>>>> 298       ierr = KSPSetOperators(ksp,J,J,SAME_NONZERO_PATTERN);
>>>> CHKERRQ(ierr);
>>>> 299       ierr = KSPSetFromOptions(ksp);    CHKERRQ(ierr);
>>>> 300       ierr = KSPGetPC(ksp,&pc);         CHKERRQ(ierr);
>>>> 301       ierr = PCFactorSetFill(pc,1);     CHKERRQ(ierr);
>>>> 302       ierr = PCFactorSetMatOrderingType(pc,"FWNOrder"); CHKERRQ(ierr);
>>>> 303       ierr = KSPSolve(ksp,b,x); CHKERRQ(ierr);
>>>> 304       ierr = KSPDestroy(&ksp);  CHKERRQ(ierr);
>>>
>>> It looks like you are putting all of this in PCApply_YourShell(). I suggest
>>> making
>>> struct YourShellContext {
>>>   PC fenpc,fwnpc;
>>>   Vec work;
>>>   PetscBool issetup;
>>> };
>>> Then pass a pointer to this thing in when you create the PCShell. Then
>>> something like
>>> PCSetUp_YourShell(PC pc) {
>>> struct YourShellContext *ctx;
>>> Mat A,B;
>>> MatStructure mstr;
>>> const char *prefix;
>>> char iprefix[256];
>>> PCShellGetContext(pc,&ctx);
>>> PCGetOptionsPrefix(pc,&prefix);
>>> if (!ctx->issetup) {
>>>   ierr = MatGetVecs(A,&ctx->work,PETSC_NULL);
>>>   PCCreate(comm,&ctx->fenpc);
>>>   PCSetType(ctx->fenpc,PCILU);
>>>   PCFactorSetFill(ctx->fenpc,1);
>>>   PCFactorSetMatOrderingType(ctx->fenpc,"FENOrder");
>>>   snprintf(iprefix,sizeof iprefix,"%sfen",prefix);
>>>   PCSetOptionsPrefix(ctx->fenpc,iprefix);
>>>   PCSetFromOptions(ctx->fenpc);
>>>   /* Same as above for ctx->fwnpc */
>>> }
>>> PCGetOperators(pc,&A,&B,&mstr);
>>> PCSetOperators(ctx->fenpc,A,B,mstr);
>>> PCSetOperators(ctx->fwnpc,A,B,mstr);
>>> PCSetUp(ctx->fenpc);
>>> PCSetUp(ctx->fwnpc);
>>> ctx->issetup = PETSC_TRUE;
>>> }
>>> Then in PCApply_YourShell(PC pc,Vec X,Vec Y) {
>>> ...
>>> PCGetOperators(pc,&A,PETSC_NULL,PETSC_NULL);
>>> PCApply(ctx->fenpc,X,Y); /* apply first ordering */
>>> VecScale(Y,-1.0);
>>> MatMultAdd(A,Y,X,ctx->work); /* Compute fresh residual */
>>> PCApply(ctx->fwnpc,ctx->work,Y);
>>> }
>>>>
>>>> I don't know if this is what you intended... is there a way to make
>>>> KSP recompute the factorization with the different ordering without
>>>> destroying and recreating the KSP?
>>>
>>> Store both as in the code above.
>>>
>>>>
>>>> Concerning Barry's option, js it possible to use the PETSc ilu
>>>> factorization function (and eventually the ones provided by external
>>>> libraries) inside my PCSHELL, passing it the various ordering I will
>>>> need? If so, what's the best way to access it?
>>>
>>> Does the code above answer this question?
>>
>>
>>
>> --
>> "[Je pense que] l'homme est un monde qui vaut des fois les mondes et
>> que les plus ardentes ambitions sont celles qui ont eu l'orgueil de
>> l'Anonymat" -- Non omnibus, sed mihi et tibi
>> Amedeo Modigliani
>
>



-- 
"[Je pense que] l'homme est un monde qui vaut des fois les mondes et
que les plus ardentes ambitions sont celles qui ont eu l'orgueil de
l'Anonymat" -- Non omnibus, sed mihi et tibi
Amedeo Modigliani