<div class="gmail_quote">On Fri, Jul 1, 2011 at 07:11, Gianluca Meneghello <span dir="ltr"><<a href="mailto:gianmail@gmail.com">gianmail@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
<div id=":1kw">I've tried to implement Jed's solution as it was looking faster to<br>
try. This is what I came up with for two different orders, FENOrder<br>
and FWNOrder respectively<br>
<br>
283 KSP ksp;<br>
284 PC pc;<br>
285<br>
286 ierr = KSPCreate(PETSC_COMM_WORLD,&ksp); CHKERRQ(ierr);<br>
287 ierr = KSPSetOptionsPrefix(ksp,"mg_"); CHKERRQ(ierr);<br>
288 ierr = KSPSetOperators(ksp,J,J,SAME_NONZERO_PATTERN); CHKERRQ(ierr);<br>
289 ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr);<br>
290 ierr = KSPGetPC(ksp,&pc); CHKERRQ(ierr);<br>
291 ierr = PCFactorSetFill(pc,1); CHKERRQ(ierr);<br>
292 ierr = PCFactorSetMatOrderingType(pc,"FENOrder"); CHKERRQ(ierr);<br>
293 ierr = KSPSolve(ksp,b,x); CHKERRQ(ierr);<br>
294 ierr = KSPDestroy(&ksp); CHKERRQ(ierr);<br>
295<br>
296 ierr = KSPCreate(PETSC_COMM_WORLD,&ksp); CHKERRQ(ierr);<br>
297 ierr = KSPSetOptionsPrefix(ksp,"mg_"); CHKERRQ(ierr);<br>
298 ierr = KSPSetOperators(ksp,J,J,SAME_NONZERO_PATTERN); CHKERRQ(ierr);<br>
299 ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr);<br>
300 ierr = KSPGetPC(ksp,&pc); CHKERRQ(ierr);<br>
301 ierr = PCFactorSetFill(pc,1); CHKERRQ(ierr);<br>
302 ierr = PCFactorSetMatOrderingType(pc,"FWNOrder"); CHKERRQ(ierr);<br>
303 ierr = KSPSolve(ksp,b,x); CHKERRQ(ierr);<br>
304 ierr = KSPDestroy(&ksp); CHKERRQ(ierr);<br></div></blockquote><div><br></div><div>It looks like you are putting all of this in PCApply_YourShell(). I suggest making </div><div><br></div><div>struct YourShellContext {</div>
<div> PC fenpc,fwnpc;</div><div> Vec work;</div><div> PetscBool issetup;</div><div>};</div><div><br></div><div>Then pass a pointer to this thing in when you create the PCShell. Then something like</div><div><br></div><div>
PCSetUp_YourShell(PC pc) {</div><div>struct YourShellContext *ctx;</div><div>Mat A,B;</div><div>MatStructure mstr;</div><div>const char *prefix;</div><div>char iprefix[256];</div><div><br></div><div>PCShellGetContext(pc,&ctx);</div>
<div>PCGetOptionsPrefix(pc,&prefix);</div><div><br></div><div>if (!ctx->issetup) {</div><div> ierr = MatGetVecs(A,&ctx->work,PETSC_NULL);</div><div> PCCreate(comm,&ctx->fenpc);</div><div> PCSetType(ctx->fenpc,PCILU);</div>
<div> PCFactorSetFill(ctx->fenpc,1);</div><div> PCFactorSetMatOrderingType(ctx->fenpc,"FENOrder");</div><div> snprintf(iprefix,sizeof iprefix,"%sfen",prefix);</div><div> PCSetOptionsPrefix(ctx->fenpc,iprefix);</div>
<div> PCSetFromOptions(ctx->fenpc);</div><div> /* Same as above for ctx->fwnpc */</div><div>}</div><div>PCGetOperators(pc,&A,&B,&mstr);</div><div><div>PCSetOperators(ctx->fenpc,A,B,mstr);</div><div>
PCSetOperators(ctx->fwnpc,A,B,mstr);</div></div><div>PCSetUp(ctx->fenpc);</div><div>PCSetUp(ctx->fwnpc);</div><div>ctx->issetup = PETSC_TRUE;</div><div>}</div><div><br></div><div>Then in PCApply_YourShell(PC pc,Vec X,Vec Y) {</div>
<div>...</div><div>PCGetOperators(pc,&A,PETSC_NULL,PETSC_NULL);</div><div>PCApply(ctx->fenpc,X,Y); /* apply first ordering */</div><div>VecScale(Y,-1.0);</div><div>MatMultAdd(A,Y,X,ctx->work); /* Compute fresh residual */</div>
<div>PCApply(ctx->fwnpc,ctx->work,Y);</div><div>}</div><div><br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;"><div id=":1kw">
<br>
I don't know if this is what you intended... is there a way to make<br>
KSP recompute the factorization with the different ordering without<br>
destroying and recreating the KSP?<br></div></blockquote><div><br></div><div>Store both as in the code above.</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
<div id=":1kw">
<br>
Concerning Barry's option, js it possible to use the PETSc ilu<br>
factorization function (and eventually the ones provided by external<br>
libraries) inside my PCSHELL, passing it the various ordering I will<br>
need? If so, what's the best way to access it?</div></blockquote></div><br><div>Does the code above answer this question?</div>