[petsc-dev] FieldSplit fixes in 3.3

[Jed Brown]

On Sun, Aug 19, 2012 at 2:31 AM, Dmitry Karpeev <karpeev at mcs.anl.gov> wrote:

> I just pushed a couple of fixes to PCFIELDSPLIT in petsc-3.3 (
> http://petsc.cs.iit.edu/petsc/releases/petsc-3.3/rev/943dde820f7f,


Yikes:

--- a/src/ksp/pc/impls/fieldsplit/fieldsplit.c
+++ b/src/ksp/pc/impls/fieldsplit/fieldsplit.c
@@ -14,6 +14,7 @@
   PetscInt          *fields,*fields_col;
   VecScatter        sctx;
   IS                is,is_col;
+  DM                dm;
   PC_FieldSplitLink next,previous;
 };

@@ -263,6 +264,13 @@
       if(dms) {
         ierr = PetscInfo(pc,"Setting up physics based fieldsplit
preconditioner using the embedded DM\n");CHKERRQ(ierr);
         for (ilink=jac->head,i=0; ilink; ilink=ilink->next,i++) {
+          /*
+           HACK:
+           keep a handle to DM here without increfing it: may need the DM
to set up the Schur solvers;
+           can't rely on KSPGetDM() since it will create a DMShell if none
was set;
+           can't use ilink->ksp->dm without creating a dependence on
dmimpl.h.
+           */
+          ilink->dm = dms[i];

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 *capability*, not *presence*. PETSc should never do
something different based on existence of a DM.

-      /* set tabbing and options prefix of KSP inside the MatSchur */
+      /* set tabbing, options prefix and DM of KSP inside the MatSchur */
       ierr  = MatSchurComplementGetKSP(jac->schur,&ksp);CHKERRQ(ierr);
       ierr  =
PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,2);CHKERRQ(ierr);
+      {
+        PC pcinner;
+        ierr           = KSPGetPC(ksp, &pcinner); CHKERRQ(ierr);
+        ierr           =
PetscObjectIncrementTabLevel((PetscObject)pcinner,(PetscObject)pc,2);CHKERRQ(ierr);
+      }

There is KSPIncrementTabLevel() for a reason.

-      ierr = PetscSNPrintf(schurprefix,sizeof
schurprefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr);
-      ierr  = KSPSetOptionsPrefix(jac->kspschur,schurprefix);CHKERRQ(ierr);
       /* really want setfromoptions called in
PCSetFromOptions_FieldSplit(), but it is not ready yet */
       /* need to call this every time, since the jac->kspschur is freshly
created, otherwise its options never get set */
-      ierr = KSPSetFromOptions(jac->kspschur);CHKERRQ(ierr);
+      ierr = PetscSNPrintf(schurprefix,sizeof
schurprefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr);
+      ierr = KSPSetOptionsPrefix(jac->kspschur,schurprefix);CHKERRQ(ierr);
+      ierr = KSPSetFromOptions(jac->kspschur); CHKERRQ(ierr);
+      ierr = KSPSetUp(jac->kspschur);          CHKERRQ(ierr);

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

-  } else {
-    /* set up the individual PCs */
-    i    = 0;
-    ilink = jac->head;
-    while (ilink) {
-      ierr =
KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i],flag);CHKERRQ(ierr);
-      /* really want setfromoptions called in
PCSetFromOptions_FieldSplit(), but it is not ready yet */
-      if (!jac->suboptionsset) {ierr =
KSPSetFromOptions(ilink->ksp);CHKERRQ(ierr);}
-      ierr = KSPSetUp(ilink->ksp);CHKERRQ(ierr);
-      i++;
-      ilink = ilink->next;
-    }
+  }
+  /* set up the individual PCs */
+  i    = 0;
+  ilink = jac->head;
+  while (ilink) {
+    ierr = KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i],flag);
CHKERRQ(ierr);
+    /* really want setfromoptions called in PCSetFromOptions_FieldSplit(),
but it is not ready yet */
+    if (!jac->suboptionsset) {ierr = KSPSetFromOptions(ilink->ksp);
CHKERRQ(ierr);}
+    ierr = KSPSetUp(ilink->ksp);CHKERRQ(ierr);
+    i++;
+    ilink = ilink->next;

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?


> http://petsc.cs.iit.edu/petsc/releases/petsc-3.3/rev/5a9ebf885615)
>

       /* HACK: Check for the constant null space */
-      ierr = MatGetNullSpace(pc->pmat, &sp);CHKERRQ(ierr);
+      /*
+       This replaces a previous test that checked whether the split's
constants are injected into the ambient problem's nullspace.
+       That can create false negatives (declaring constants not in the
nullspace, when they are) as in Stokes or A = [1, 1; 1, 0].
+       */
       if (sp) {

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.

       ierr  = MatGetNullSpace(jac->pmat[1], &sp);CHKERRQ(ierr);
-      if (sp) {ierr  = MatSetNullSpace(jac->schur, sp);CHKERRQ(ierr);}
-      /* set tabbing, options prefix and DM of KSP inside the MatSchur */
+      /*
+       Do we really want to attach the A11-block's nullspace to S? Note
that this may generate
+       "false positives", when the A11's nullspace isn't S's: Stokes or A
= [1, 1; 1, 0].
+       Using a false nullspace may prevent KSP(S) from converging, since
 it might force an inconsistent rhs.
+       */
+      /* if (sp) {ierr  = MatSetNullSpace(jac->schur, sp);CHKERRQ(ierr);}
*/
+      /* set tabbing, options prefix and DM of KSP inside the
MatSchurComplement */

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.


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.

 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.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20120819/16ea96d7/attachment.html>