[petsc-users] ILU and Block Gauss Seidel smoothing

[Barry Smith]

On Jul 13, 2011, at 4:02 AM, Gianluca Meneghello wrote:

> 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.

   Definitely not. Looking at PCSetUp_ILU() 

 } else {
    if (!pc->setupcalled) {
      /* first time in so compute reordering and symbolic factorization */
      ierr = MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);CHKERRQ(ierr);

it will actually call the ordering routine a TOTAL of ONE TIME if you use SAME_NONZERO_PATTERN to KSPSetOperators() or PCSetOperators(). (it is smart enough to reuse the ordering).

   BTW: I don't think you even need to use a PCSHELL preconditioner. You can get the same effect with a PCCOMPOSITE.  Just create a single KSP in the usual way and use the options

-pc_type composite -pc_composite_type ilu,ilu -sub_0_ksp_type preonly -sub_1_ksp_type preonly -sub_0_pc_factor_mat_ordering_type FESOrder -sub_1_pc_factor_mat_ordering_type anotherorder somethingelse


   Barry





> 
> 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