<div class="gmail_quote">On Sun, Aug 19, 2012 at 2:31 AM, Dmitry Karpeev <span dir="ltr"><<a href="mailto:karpeev@mcs.anl.gov" target="_blank">karpeev@mcs.anl.gov</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
I just pushed a couple of fixes to PCFIELDSPLIT in petsc-3.3 (<a href="http://petsc.cs.iit.edu/petsc/releases/petsc-3.3/rev/943dde820f7f" target="_blank">http://petsc.cs.iit.edu/petsc/releases/petsc-3.3/rev/943dde820f7f</a>, </blockquote>
<div><br></div><div>Yikes:</div><div><br></div><div><div><font face="courier new, monospace">--- a/src/ksp/pc/impls/fieldsplit/fieldsplit.c</font></div><div><font face="courier new, monospace">+++ b/src/ksp/pc/impls/fieldsplit/fieldsplit.c</font></div>
<div><font face="courier new, monospace">@@ -14,6 +14,7 @@</font></div><div><font face="courier new, monospace"> PetscInt *fields,*fields_col;</font></div><div><font face="courier new, monospace"> VecScatter sctx;</font></div>
<div><font face="courier new, monospace"> IS is,is_col;</font></div><div><font face="courier new, monospace">+ DM dm;</font></div><div><font face="courier new, monospace"> PC_FieldSplitLink next,previous;</font></div>
<div><font face="courier new, monospace"> };</font></div><div><font face="courier new, monospace"> </font></div><div><font face="courier new, monospace">@@ -263,6 +264,13 @@</font></div><div><font face="courier new, monospace"> if(dms) {</font></div>
<div><font face="courier new, monospace"> ierr = PetscInfo(pc,"Setting up physics based fieldsplit preconditioner using the embedded DM\n");CHKERRQ(ierr);</font></div><div><font face="courier new, monospace"> for (ilink=jac->head,i=0; ilink; ilink=ilink->next,i++) {</font></div>
<div><font face="courier new, monospace">+ /* </font></div><div><font face="courier new, monospace">+ HACK: </font></div><div><font face="courier new, monospace">+ keep a handle to DM here without increfing it: may need the DM to set up the Schur solvers; </font></div>
<div><font face="courier new, monospace">+ can't rely on KSPGetDM() since it will create a DMShell if none was set;</font></div><div><font face="courier new, monospace">+ can't use ilink->ksp->dm without creating a dependence on dmimpl.h.</font></div>
<div><font face="courier new, monospace">+ */</font></div><div><font face="courier new, monospace">+ ilink->dm = dms[i]; </font></div></div><div><br></div><div>I think your comment was intended to say "kspimpl.h", but you are creating a tangled web of barbed wire here. You are supposed to call KSPGetDM() and switch based on its <b>capability</b>, not <b>presence</b>. PETSc should never do something different based on existence of a DM.</div>
<div><br></div><div><div><font face="courier new, monospace">- /* set tabbing and options prefix of KSP inside the MatSchur */</font></div><div><font face="courier new, monospace">+ /* set tabbing, options prefix and DM of KSP inside the MatSchur */</font></div>
<div><font face="courier new, monospace"> ierr = MatSchurComplementGetKSP(jac->schur,&ksp);CHKERRQ(ierr);</font></div><div><font face="courier new, monospace"> ierr = PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,2);CHKERRQ(ierr);</font></div>
<div><font face="courier new, monospace">+ {</font></div><div><font face="courier new, monospace">+ PC pcinner;</font></div><div><font face="courier new, monospace">+ ierr = KSPGetPC(ksp, &pcinner); CHKERRQ(ierr);</font></div>
<div><font face="courier new, monospace">+ ierr = PetscObjectIncrementTabLevel((PetscObject)pcinner,(PetscObject)pc,2);CHKERRQ(ierr);</font></div><div><font face="courier new, monospace">+ }</font></div>
</div><div><br></div><div>There is KSPIncrementTabLevel() for a reason.</div><div><br></div><div><div><font face="courier new, monospace">- ierr = PetscSNPrintf(schurprefix,sizeof schurprefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr);</font></div>
<div><font face="courier new, monospace">- ierr = KSPSetOptionsPrefix(jac->kspschur,schurprefix);CHKERRQ(ierr);</font></div><div><font face="courier new, monospace"> /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */</font></div>
<div><font face="courier new, monospace"> /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */</font></div><div><font face="courier new, monospace">- ierr = KSPSetFromOptions(jac->kspschur);CHKERRQ(ierr);</font></div>
<div><font face="courier new, monospace">+ ierr = PetscSNPrintf(schurprefix,sizeof schurprefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr);</font></div>
<div><font face="courier new, monospace">+ ierr = KSPSetOptionsPrefix(jac->kspschur,schurprefix);CHKERRQ(ierr);</font></div><div><font face="courier new, monospace">+ ierr = KSPSetFromOptions(jac->kspschur); CHKERRQ(ierr);</font></div>
<div><font face="courier new, monospace">+ ierr = KSPSetUp(jac->kspschur); CHKERRQ(ierr);</font></div></div><div><br></div><div>Why this KSPSetUp()? You have changed the formatting here, but not explained why. (Also, can't you follow PETSc coding guidelines regarding placement of CHKERRQ?)</div>
<div><br></div><div><div><font face="courier new, monospace">- } else {</font></div><div><font face="courier new, monospace">- /* set up the individual PCs */</font></div><div><font face="courier new, monospace">- i = 0;</font></div>
<div><font face="courier new, monospace">- ilink = jac->head;</font></div><div><font face="courier new, monospace">- while (ilink) {</font></div><div><font face="courier new, monospace">- ierr = KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i],flag);CHKERRQ(ierr);</font></div>
<div><font face="courier new, monospace">- /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */</font></div><div><font face="courier new, monospace">- if (!jac->suboptionsset) {ierr = KSPSetFromOptions(ilink->ksp);CHKERRQ(ierr);}</font></div>
<div><font face="courier new, monospace">- ierr = KSPSetUp(ilink->ksp);CHKERRQ(ierr);</font></div><div><font face="courier new, monospace">- i++;</font></div><div><font face="courier new, monospace">- ilink = ilink->next;</font></div>
<div><font face="courier new, monospace">- }</font></div><div><font face="courier new, monospace">+ } </font></div><div><font face="courier new, monospace">+ /* set up the individual PCs */</font></div><div><font face="courier new, monospace">+ i = 0;</font></div>
<div><font face="courier new, monospace">+ ilink = jac->head;</font></div><div><font face="courier new, monospace">+ while (ilink) {</font></div><div><font face="courier new, monospace">+ ierr = KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i],flag); CHKERRQ(ierr);</font></div>
<div><font face="courier new, monospace">+ /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */</font></div><div><font face="courier new, monospace">+ if (!jac->suboptionsset) {ierr = KSPSetFromOptions(ilink->ksp); CHKERRQ(ierr);}</font></div>
<div><font face="courier new, monospace">+ ierr = KSPSetUp(ilink->ksp);CHKERRQ(ierr);</font></div><div><font face="courier new, monospace">+ i++;</font></div><div><font face="courier new, monospace">+ ilink = ilink->next;</font></div>
</div><div><br></div><div>Why do you want this stuff done even for Schur? At least one of these blocks is not used at all, so why are we setting it up?</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div><a href="http://petsc.cs.iit.edu/petsc/releases/petsc-3.3/rev/5a9ebf885615" target="_blank">http://petsc.cs.iit.edu/petsc/releases/petsc-3.3/rev/5a9ebf885615</a>)</div></blockquote><div><br></div><div><font face="courier new, monospace"> /* HACK: Check for the constant null space */</font></div>
<div><font face="courier new, monospace">- ierr = MatGetNullSpace(pc->pmat, &sp);CHKERRQ(ierr);</font></div><div><div><font face="courier new, monospace">+ /* </font></div><div><font face="courier new, monospace">+ This replaces a previous test that checked whether the split's constants are injected into the ambient problem's nullspace.</font></div>
<div><font face="courier new, monospace">+ That can create false negatives (declaring constants not in the nullspace, when they are) as in Stokes or A = [1, 1; 1, 0].</font></div><div><font face="courier new, monospace">+ */</font></div>
<div><font face="courier new, monospace"> if (sp) {</font></div></div><div><br></div><div>So now you have stale comments and conditionals based on uninitialized memory (sp)? I suppose your intent is to always do this constant null space check regardless of whether the larger matrix has a null space set. I agree that the old code was incorrect because the restriction of the null space of the parent matrix is not comparable to the null space of the submatrix.</div>
<div><br></div><div><div><font face="courier new, monospace"> ierr = MatGetNullSpace(jac->pmat[1], &sp);CHKERRQ(ierr);</font></div><div><font face="courier new, monospace">- if (sp) {ierr = MatSetNullSpace(jac->schur, sp);CHKERRQ(ierr);}</font></div>
<div><font face="courier new, monospace">- /* set tabbing, options prefix and DM of KSP inside the MatSchur */</font></div><div><font face="courier new, monospace">+ /* </font></div><div><font face="courier new, monospace">+ Do we really want to attach the A11-block's nullspace to S? Note that this may generate</font></div>
<div><font face="courier new, monospace">+ "false positives", when the A11's nullspace isn't S's: Stokes or A = [1, 1; 1, 0].</font></div><div><font face="courier new, monospace">+ Using a false nullspace may prevent KSP(S) from converging, since it might force an inconsistent rhs. </font></div>
<div><font face="courier new, monospace">+ */</font></div><div><font face="courier new, monospace">+ /* if (sp) {ierr = MatSetNullSpace(jac->schur, sp);CHKERRQ(ierr);} */</font></div><div><font face="courier new, monospace">+ /* set tabbing, options prefix and DM of KSP inside the MatSchurComplement */</font></div>
</div><div><br></div><div>I agree with this change, though the comment is incorrect. The typical problem is that S is nonsingular or has a small null space while A11 is has huge null space (e.g. the zero matrix). The true null space projection for the zero matrix prevents progress in any direction.</div>
<div><br></div><div><br></div><div>I think we can repair the second patch, though we really need a general mechanism for providing Schur complement null spaces. Matt advocates attaching it to A11 (where the user knows full well that it's not the null space of A11 when they set it), making the old code "correct". I think that's confusing and we need a specific API for it, but if people are using it that way, we shouldn't change it in 3.3.</div>
<div><br></div><div> I'd prefer to backout the first patch, split it into separate pieces, and review each before pushing. I think it's a potentially big behavior change for 3.3.</div></div>