Barry,
thanks again for your input. That would be really good.
When executing my program with the options you suggested I obtain
various errors (unknown pc etc) so I tried to figure out how to avoid
them. Is it something like this what you intended?
./main
-mg_pc_type composite
-mg_pc_composite_pcs ilu,ilu
-mg_ksp_type preonly
-mg_sub_0_pc_factor_mat_ordering_type FESOrder
-mg_sub_1_pc_factor_mat_ordering_type FWSOrder
-mg_pc_composite_type multiplicative
-mg_ksp_view
provided that
KSPSetOptionsPrefix(ksp,"mg_");
has been applied to the KSP I want to modify.
Thanks again
Gianluca
On 13 July 2011 15:33, Barry Smith <bsmith at mcs.anl.gov> wrote:
>
> 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
>
>
--
"[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
