Commits

Matt Knepley committed 2e71c61 Merge with conflicts

Merge remote-tracking branch 'origin/karpeev/ksp-matschurcomplement'

* origin/karpeev/ksp-matschurcomplement:
Fixed memory leak in MatSchurComplement.
MatSchurComplementXXXSubMatrices renamed for consistency.
PCFieldSplit: enable extraction of MatSchurComplement to configure it.
MatSchurComplement can use different approximations to inv(A00) when forming Sp.
PCFieldSplit: Schur can use lumping with MatSchurComplement.
MatSchurComplement: Allow lumping of A00 when forming Sp.
PCFieldSplit implements SELFP preconditioning for Schur using MatSchurComplementGetPmat to assemble Sp.
MatCreateSchurComplementPmat() now public.
MatSchurComplement now uses MatDiagonalScale to form Sp. Factored out Sp constructor for consistency and reuse.
More consistent matrix naming in arguments and docs to MatSchurComplement-related routines.
MatSchurComplement can now assemble Sp on demand.

Conflicts:
include/petscksp.h
include/petscpc.h
src/ksp/pc/impls/fieldsplit/fieldsplit.c

  • Participants
  • Parent commits 1a77d57, 3b1e2a3

Comments (0)

Files changed (7)

File include/petscksp.h

 PETSC_EXTERN PetscErrorCode KSPInitializePackage(void);
 
 /*S
-     KSP - Abstract PETSc object that manages all Krylov methods. This is the object that manages the 
+     KSP - Abstract PETSc object that manages all Krylov methods. This is the object that manages the
          linear solves in PETSc (even those such as direct solvers that do no use Krylov accelerators).
 
    Level: beginner
 PETSC_EXTERN PetscErrorCode KSPSetFischerGuess(KSP,KSPFischerGuess);
 PETSC_EXTERN PetscErrorCode KSPGetFischerGuess(KSP,KSPFischerGuess*);
 
+/*E
+    MatSchurComplementAinvType - Determines how to approximate the inverse of the (0,0) block in Schur complement preconditioning matrix assembly routines
+
+    Level: intermediate
+
+.seealso: MatSchurComplementGetAinvType(), MatSchurComplementSetAinvType(), MatSchurComplementGetPmat(), MatGetSchurComplement(), MatCreateSchurComplementPmat()
+E*/
+typedef enum {MAT_SCHUR_COMPLEMENT_AINV_DIAG, MAT_SCHUR_COMPLEMENT_AINV_LUMP} MatSchurComplementAinvType;
+PETSC_EXTERN const char *const MatSchurComplementAinvTypes[];
+
 PETSC_EXTERN PetscErrorCode MatCreateSchurComplement(Mat,Mat,Mat,Mat,Mat,Mat*);
 PETSC_EXTERN PetscErrorCode MatSchurComplementGetKSP(Mat,KSP*);
 PETSC_EXTERN PetscErrorCode MatSchurComplementSetKSP(Mat,KSP);
-PETSC_EXTERN PetscErrorCode MatSchurComplementSet(Mat,Mat,Mat,Mat,Mat,Mat);
-PETSC_EXTERN PetscErrorCode MatSchurComplementUpdate(Mat,Mat,Mat,Mat,Mat,Mat,MatStructure);
-PETSC_EXTERN PetscErrorCode MatSchurComplementGetSubmatrices(Mat,Mat*,Mat*,Mat*,Mat*,Mat*);
+PETSC_EXTERN PetscErrorCode MatSchurComplementSetSubMatrices(Mat,Mat,Mat,Mat,Mat,Mat);
+PETSC_EXTERN PetscErrorCode MatSchurComplementUpdateSubMatrices(Mat,Mat,Mat,Mat,Mat,Mat,MatStructure);
+PETSC_EXTERN PetscErrorCode MatSchurComplementGetSubMatrices(Mat,Mat*,Mat*,Mat*,Mat*,Mat*);
+PETSC_EXTERN PetscErrorCode MatSchurComplementSetAinvType(Mat,MatSchurComplementAinvType);
+PETSC_EXTERN PetscErrorCode MatSchurComplementGetAinvType(Mat,MatSchurComplementAinvType*);
+PETSC_EXTERN PetscErrorCode MatSchurComplementGetPmat(Mat,MatReuse,Mat*);
 PETSC_EXTERN PetscErrorCode MatSchurComplementComputeExplicitOperator(Mat,Mat*);
-PETSC_EXTERN PetscErrorCode MatGetSchurComplement(Mat,IS,IS,IS,IS,MatReuse,Mat *,MatReuse,Mat *);
+PETSC_EXTERN PetscErrorCode MatGetSchurComplement(Mat,IS,IS,IS,IS,MatReuse,Mat *,MatSchurComplementAinvType,MatReuse,Mat *);
+PETSC_EXTERN PetscErrorCode MatCreateSchurComplementPmat(Mat,Mat,Mat,Mat,MatSchurComplementAinvType,MatReuse,Mat*);
 
 PETSC_EXTERN PetscErrorCode KSPSetDM(KSP,DM);
 PETSC_EXTERN PetscErrorCode KSPSetDMActive(KSP,PetscBool );

File include/petscpc.h

 
 .seealso: PCFieldSplitSchurPrecondition()
 E*/
-typedef enum {PC_FIELDSPLIT_SCHUR_PRE_SELF,PC_FIELDSPLIT_SCHUR_PRE_A11,PC_FIELDSPLIT_SCHUR_PRE_USER,PC_FIELDSPLIT_SCHUR_PRE_FULL} PCFieldSplitSchurPreType;
+typedef enum {PC_FIELDSPLIT_SCHUR_PRE_SELF,PC_FIELDSPLIT_SCHUR_PRE_SELFP,PC_FIELDSPLIT_SCHUR_PRE_A11,PC_FIELDSPLIT_SCHUR_PRE_USER,PC_FIELDSPLIT_SCHUR_PRE_FULL} PCFieldSplitSchurPreType;
 PETSC_EXTERN const char *const PCFieldSplitSchurPreTypes[];
 
 /*E
 PETSC_EXTERN PetscErrorCode PCFieldSplitSchurPrecondition(PC,PCFieldSplitSchurPreType,Mat);
 PETSC_EXTERN PetscErrorCode PCFieldSplitSetSchurFactType(PC,PCFieldSplitSchurFactType);
 PETSC_EXTERN PetscErrorCode PCFieldSplitGetSchurBlocks(PC,Mat*,Mat*,Mat*,Mat*);
+PETSC_EXTERN PetscErrorCode PCFieldSplitSchurGetS(PC,Mat *S);
+PETSC_EXTERN PetscErrorCode PCFieldSplitSchurRestoreS(PC,Mat *S);
 
 PETSC_EXTERN PetscErrorCode PCGalerkinSetRestriction(PC,Mat);
 PETSC_EXTERN PetscErrorCode PCGalerkinSetInterpolation(PC,Mat);

File src/ksp/ksp/examples/tests/ex21.c

   ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
   ierr = ISView(is0,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
   ierr = ISView(is1,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
-  ierr = MatGetSchurComplement(A,is0,is0,is1,is1,MAT_INITIAL_MATRIX,&S,MAT_IGNORE_MATRIX,NULL);CHKERRQ(ierr);
+  ierr = MatGetSchurComplement(A,is0,is0,is1,is1,MAT_INITIAL_MATRIX,&S,PETSC_FALSE,MAT_IGNORE_MATRIX,NULL);CHKERRQ(ierr);
   ierr = MatComputeExplicitOperator(S,&Sexplicit);CHKERRQ(ierr);
   ierr = PetscPrintf(PETSC_COMM_WORLD,"\nExplicit Schur complement of (0,0) in (1,1)\n");CHKERRQ(ierr);
   ierr = MatView(Sexplicit,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
 
   /* And the other */
   ierr = Create(PETSC_COMM_WORLD,&A,&is0,&is1);CHKERRQ(ierr);
-  ierr = MatGetSchurComplement(A,is1,is1,is0,is0,MAT_INITIAL_MATRIX,&S,MAT_IGNORE_MATRIX,NULL);CHKERRQ(ierr);
+  ierr = MatGetSchurComplement(A,is1,is1,is0,is0,MAT_INITIAL_MATRIX,&S,PETSC_FALSE,MAT_IGNORE_MATRIX,NULL);CHKERRQ(ierr);
   ierr = MatComputeExplicitOperator(S,&Sexplicit);CHKERRQ(ierr);
   ierr = PetscPrintf(PETSC_COMM_WORLD,"\nExplicit Schur complement of (1,1) in (0,0)\n");CHKERRQ(ierr);
   ierr = MatView(Sexplicit,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
 
   /* This time just the preconditioner */
   ierr = Create(PETSC_COMM_WORLD,&A,&is0,&is1);CHKERRQ(ierr);
-  ierr = MatGetSchurComplement(A,is0,is0,is1,is1,MAT_IGNORE_MATRIX,NULL,MAT_INITIAL_MATRIX,&S);CHKERRQ(ierr);
+  ierr = MatGetSchurComplement(A,is0,is0,is1,is1,MAT_IGNORE_MATRIX,NULL,PETSC_FALSE,MAT_INITIAL_MATRIX,&S);CHKERRQ(ierr);
   ierr = PetscPrintf(PETSC_COMM_WORLD,"\nPreconditioning Schur complement of (0,0) in (1,1)\n");CHKERRQ(ierr);
   ierr = MatView(S,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
   /* Modify and refresh */
   ierr = MatShift(A,1.);CHKERRQ(ierr);
-  ierr = MatGetSchurComplement(A,is0,is0,is1,is1,MAT_IGNORE_MATRIX,NULL,MAT_REUSE_MATRIX,&S);CHKERRQ(ierr);
+  ierr = MatGetSchurComplement(A,is0,is0,is1,is1,MAT_IGNORE_MATRIX,NULL,PETSC_FALSE,MAT_REUSE_MATRIX,&S);CHKERRQ(ierr);
   ierr = PetscPrintf(PETSC_COMM_WORLD,"\nAfter update\n");CHKERRQ(ierr);
   ierr = MatView(S,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
   ierr = Destroy(&A,&is0,&is1);CHKERRQ(ierr);

File src/ksp/ksp/utils/schurm.c

 
 #include <petsc-private/matimpl.h>
 #include <petscksp.h>                 /*I "petscksp.h" I*/
+const char *const MatSchurComplementAinvTypes[] = {"SELF","LUMP","MatSchurComplementAinvType","MAT_SCHUR_COMPLEMENT_AINV_",0};
 
 typedef struct {
-  Mat A,Ap,B,C,D;
-  KSP ksp;
-  Vec work1,work2;
+  Mat                        A,Ap,B,C,D;
+  KSP                        ksp;
+  Vec                        work1,work2;
+  MatSchurComplementAinvType ainvtype;
 } Mat_SchurComplement;
 
+
+
 #undef __FUNCT__
 #define __FUNCT__ "MatGetVecs_SchurComplement"
 PetscErrorCode MatGetVecs_SchurComplement(Mat N,Vec *right,Vec *left)
 
 
 /*
-           A11 - A10 ksp(A00,Ap00)  A01
+           A11 - A10 ksp(A00,Ap00) A01
 */
 #undef __FUNCT__
 #define __FUNCT__ "MatMult_SchurComplement"
 }
 
 /*
-           A11 - A10 ksp(A00,Ap00)  A01
+           A11 - A10 ksp(A00,Ap00) A01
 */
 #undef __FUNCT__
 #define __FUNCT__ "MatMultAdd_SchurComplement"
   PetscErrorCode      ierr;
 
   PetscFunctionBegin;
+  ierr = PetscOptionsHead("MatSchurComplementOptions");CHKERRQ(ierr);
+  Na->ainvtype = MAT_SCHUR_COMPLEMENT_AINV_DIAG;
+  ierr = PetscOptionsEnum("-mat_schur_complement_ainv_type","Type of approximation for inv(A00) used when assembling Sp = A11 - A10 inv(A00) A01","MatSchurComplementSetAinvType",MatSchurComplementAinvTypes,(PetscEnum)Na->ainvtype,(PetscEnum*)&Na->ainvtype,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsTail();CHKERRQ(ierr);
   ierr = KSPSetFromOptions(Na->ksp);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
    Collective on Mat
 
-   Input Parameter:
-.   A00,A01,A10,A11  - the four parts of the original matrix (A00 is optional)
+   Input Parameters:
++   A00,A01,A10,A11  - the four parts of the original matrix A = [A00 A01; A10 A11] (A11 is optional)
+-   Ap00             - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A^{-1}
 
    Output Parameter:
-.   N - the matrix that the Schur complement A11 - A10 ksp(A00) A01
+.   S - the matrix that the Schur complement S = A11 - A10 ksp(A00,Ap00) A01
 
    Level: intermediate
 
-   Notes: The Schur complement is NOT actually formed! Rather this
-          object performs the matrix-vector product by using the the formula for
-          the Schur complement and a KSP solver to approximate the action of inv(A)
+   Notes: The Schur complement is NOT actually formed! Rather, this
+          object performs the matrix-vector product by using formula S = A11 - A10 A^{-1} A01
+	  for Schur complement S and a KSP solver to approximate the action of A^{-1}.
 
-          All four matrices must have the same MPI communicator
+          All four matrices must have the same MPI communicator.
 
-          A00 and  A11 must be square matrices
+          A00 and  A11 must be square matrices.
 
-.seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatSchurComplementUpdate(), MatCreateTranspose(), MatGetSchurComplement()
+.seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatSchurComplementUpdateSubMatrices(), MatCreateTranspose(), MatGetSchurComplement()
 
 @*/
-PetscErrorCode  MatCreateSchurComplement(Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11,Mat *N)
+PetscErrorCode  MatCreateSchurComplement(Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11,Mat *S)
 {
   PetscErrorCode ierr;
 
   PetscFunctionBegin;
   ierr = KSPInitializePackage();CHKERRQ(ierr);
-  ierr = MatCreate(((PetscObject)A00)->comm,N);CHKERRQ(ierr);
-  ierr = MatSetType(*N,MATSCHURCOMPLEMENT);CHKERRQ(ierr);
-  ierr = MatSchurComplementSet(*N,A00,Ap00,A01,A10,A11);CHKERRQ(ierr);
+  ierr = MatCreate(((PetscObject)A00)->comm,S);CHKERRQ(ierr);
+  ierr = MatSetType(*S,MATSCHURCOMPLEMENT);CHKERRQ(ierr);
+  ierr = MatSchurComplementSetSubMatrices(*S,A00,Ap00,A01,A10,A11);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
 #undef __FUNCT__
-#define __FUNCT__ "MatSchurComplementSet"
+#define __FUNCT__ "MatSchurComplementSetSubMatrices"
 /*@
-      MatSchurComplementSet - Sets the matrices that define the Schur complement
+      MatSchurComplementSetSubMatrices - Sets the matrices that define the Schur complement
 
    Collective on Mat
 
    Input Parameter:
-+   N - matrix obtained with MatCreate() and MatSetType(MATSCHURCOMPLEMENT);
--   A00,A01,A10,A11  - the four parts of the original matrix (A00 is optional)
++   S                - matrix obtained with MatCreateSchurComplement (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
+.   A00,A01,A10,A11  - the four parts of A = [A00 A01; A10 A11] (A11 is optional)
+-   Ap00             - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A^{-1}.
 
    Level: intermediate
 
-   Notes: The Schur complement is NOT actually formed! Rather this
-          object performs the matrix-vector product by using the the formula for
-          the Schur complement and a KSP solver to approximate the action of inv(A)
+   Notes: The Schur complement is NOT actually formed! Rather, this
+          object performs the matrix-vector product by using formula S = A11 - A10 A^{-1} A01
+	  for Schur complement S and a KSP solver to approximate the action of A^{-1}.
 
-          All four matrices must have the same MPI communicator
+          All four matrices must have the same MPI communicator.
 
-          A00 and  A11 must be square matrices
+          A00 and  A11 must be square matrices.
 
-.seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatSchurComplementUpdate(), MatCreateTranspose(), MatGetSchurComplement()
+.seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatSchurComplementUpdateSubMatrices(), MatCreateTranspose(), MatCreateSchurComplement(), MatGetSchurComplement()
 
 @*/
-PetscErrorCode  MatSchurComplementSet(Mat N,Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11)
+PetscErrorCode  MatSchurComplementSetSubMatrices(Mat S,Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11)
 {
   PetscErrorCode      ierr;
   PetscInt            m,n;
-  Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
+  Mat_SchurComplement *Na = (Mat_SchurComplement*)S->data;
 
   PetscFunctionBegin;
-  if (N->assembled) SETERRQ(PetscObjectComm((PetscObject)N),PETSC_ERR_ARG_WRONGSTATE,"Use MatSchurComplementUpdate() for already used matrix");
+  if (S->assembled) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Use MatSchurComplementUpdateSubMatrices() for already used matrix");
   PetscValidHeaderSpecific(A00,MAT_CLASSID,1);
   PetscValidHeaderSpecific(Ap00,MAT_CLASSID,2);
   PetscValidHeaderSpecific(A01,MAT_CLASSID,3);
 
   ierr   = MatGetLocalSize(A01,NULL,&n);CHKERRQ(ierr);
   ierr   = MatGetLocalSize(A10,&m,NULL);CHKERRQ(ierr);
-  ierr   = MatSetSizes(N,m,n,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
+  ierr   = MatSetSizes(S,m,n,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
   ierr   = PetscObjectReference((PetscObject)A00);CHKERRQ(ierr);
   ierr   = PetscObjectReference((PetscObject)Ap00);CHKERRQ(ierr);
   ierr   = PetscObjectReference((PetscObject)A01);CHKERRQ(ierr);
   if (A11) {
     ierr = PetscObjectReference((PetscObject)A11);CHKERRQ(ierr);
   }
-  N->assembled    = PETSC_TRUE;
-  N->preallocated = PETSC_TRUE;
+  S->assembled    = PETSC_TRUE;
+  S->preallocated = PETSC_TRUE;
 
-  ierr = PetscLayoutSetUp((N)->rmap);CHKERRQ(ierr);
-  ierr = PetscLayoutSetUp((N)->cmap);CHKERRQ(ierr);
+  ierr = PetscLayoutSetUp((S)->rmap);CHKERRQ(ierr);
+  ierr = PetscLayoutSetUp((S)->cmap);CHKERRQ(ierr);
   ierr = KSPSetOperators(Na->ksp,A00,Ap00,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 #undef __FUNCT__
 #define __FUNCT__ "MatSchurComplementGetKSP"
 /*@
-  MatSchurComplementGetKSP - Gets the KSP object that is used to invert A in the Schur complement matrix S = C A^{-1} B
+  MatSchurComplementGetKSP - Gets the KSP object that is used to invert A00 in the Schur complement matrix S = A11 - A10 ksp(A00,Ap00) A01
 
   Not Collective
 
   Input Parameter:
-. A - matrix created with MatCreateSchurComplement()
+. S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
 
   Output Parameter:
 . ksp - the linear solver object
 
   Options Database:
-. -fieldsplit_0_XXX sets KSP and PC options for the A block solver inside the Schur complement
+. -fieldsplit_<splitname_0>_XXX sets KSP and PC options for the 0-split solver inside the Schur complement used in PCFieldSplit; default <splitname_0> is 0.
 
   Level: intermediate
 
 .seealso: MatSchurComplementSetKSP(), MatCreateSchurComplement(), MatCreateNormal(), MatMult(), MatCreate()
 @*/
-PetscErrorCode MatSchurComplementGetKSP(Mat A, KSP *ksp)
+PetscErrorCode MatSchurComplementGetKSP(Mat S, KSP *ksp)
 {
   Mat_SchurComplement *Na;
 
   PetscFunctionBegin;
-  PetscValidHeaderSpecific(A,MAT_CLASSID,1);
+  PetscValidHeaderSpecific(S,MAT_CLASSID,1);
   PetscValidPointer(ksp,2);
-  Na   = (Mat_SchurComplement*) A->data;
+  Na   = (Mat_SchurComplement*) S->data;
   *ksp = Na->ksp;
   PetscFunctionReturn(0);
 }
 #undef __FUNCT__
 #define __FUNCT__ "MatSchurComplementSetKSP"
 /*@
-  MatSchurComplementSetKSP - Sets the KSP object that is used to invert A in the Schur complement matrix S = C A^{-1} B
+  MatSchurComplementSetKSP - Sets the KSP object that is used to invert A00 in the Schur complement matrix S = A11 - A10 ksp(A00,Ap00) A01
 
   Not Collective
 
   Input Parameters:
-+ A - matrix created with MatCreateSchurComplement()
++ S   - matrix created with MatCreateSchurComplement()
 - ksp - the linear solver object
 
   Level: developer
 
   Developer Notes:
-    There is likely no use case for this function.
+    This is used in PCFieldSplit to reuse the 0-split KSP to implement ksp(A00,Ap00) in S.
 
 .seealso: MatSchurComplementGetKSP(), MatCreateSchurComplement(), MatCreateNormal(), MatMult(), MatCreate(), MATSCHURCOMPLEMENT
 @*/
-PetscErrorCode MatSchurComplementSetKSP(Mat A, KSP ksp)
+PetscErrorCode MatSchurComplementSetKSP(Mat S, KSP ksp)
 {
   Mat_SchurComplement *Na;
   PetscErrorCode      ierr;
 
   PetscFunctionBegin;
-  PetscValidHeaderSpecific(A,MAT_CLASSID,1);
+  PetscValidHeaderSpecific(S,MAT_CLASSID,1);
   PetscValidHeaderSpecific(ksp,KSP_CLASSID,2);
-  Na      = (Mat_SchurComplement*) A->data;
+  Na      = (Mat_SchurComplement*) S->data;
   ierr    = PetscObjectReference((PetscObject)ksp);CHKERRQ(ierr);
   ierr    = KSPDestroy(&Na->ksp);CHKERRQ(ierr);
   Na->ksp = ksp;
 }
 
 #undef __FUNCT__
-#define __FUNCT__ "MatSchurComplementUpdate"
+#define __FUNCT__ "MatSchurComplementUpdateSubMatrices"
 /*@
-      MatSchurComplementUpdate - Updates the Schur complement matrix object with new submatrices
+      MatSchurComplementUpdateSubMatrices - Updates the Schur complement matrix object with new submatrices
 
    Collective on Mat
 
    Input Parameters:
-+   N - the matrix obtained with MatCreateSchurComplement()
-.   A,B,C,D  - the four parts of the original matrix (D is optional)
--   str - either SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER
++   S                - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
+.   A00,A01,A10,A11  - the four parts of A = [A00 A01; A10 A11] (A11 is optional)
+.   Ap00             - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A^{-1}.
+-   str              - either SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER
 
 
    Level: intermediate
 
    Notes: All four matrices must have the same MPI communicator
 
-          A and  D must be square matrices
+          A00 and  A11 must be square matrices
 
-          All of the matrices provided must have the same sizes as was used with MatCreateSchurComplement() or MatSchurComplementSet()
-          though they need not be the same matrices
+          All of the matrices provided must have the same sizes as was used with MatCreateSchurComplement() or MatSchurComplementSetSubMatrices()
+          though they need not be the same matrices.
 
 .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatCreateSchurComplement()
 
 @*/
-PetscErrorCode  MatSchurComplementUpdate(Mat N,Mat A,Mat Ap,Mat B,Mat C,Mat D,MatStructure str)
+PetscErrorCode  MatSchurComplementUpdateSubMatrices(Mat S,Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11,MatStructure str)
 {
   PetscErrorCode      ierr;
-  Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
+  Mat_SchurComplement *Na = (Mat_SchurComplement*)S->data;
 
   PetscFunctionBegin;
-  if (!N->assembled) SETERRQ(PetscObjectComm((PetscObject)N),PETSC_ERR_ARG_WRONGSTATE,"Use MatSchurComplementSet() for new matrix");
-  PetscValidHeaderSpecific(A,MAT_CLASSID,1);
-  PetscValidHeaderSpecific(B,MAT_CLASSID,2);
-  PetscValidHeaderSpecific(C,MAT_CLASSID,3);
-  PetscCheckSameComm(A,1,Ap,2);
-  PetscCheckSameComm(A,1,B,3);
-  PetscCheckSameComm(A,1,C,4);
-  if (A->rmap->n != A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of A %D do not equal local columns %D",A->rmap->n,A->cmap->n);
-  if (A->rmap->n != Ap->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of A %D do not equal local rows of Ap %D",A->rmap->n,Ap->rmap->n);
-  if (Ap->rmap->n != Ap->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of Ap %D do not equal local columns %D",Ap->rmap->n,Ap->cmap->n);
-  if (A->cmap->n != B->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local columns of A %D do not equal local rows of B %D",A->cmap->n,B->rmap->n);
-  if (C->cmap->n != A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local columns of C %D do not equal local rows of A %D",C->cmap->n,A->rmap->n);
-  if (D) {
-    PetscValidHeaderSpecific(D,MAT_CLASSID,5);
-    PetscCheckSameComm(A,1,D,5);
-    if (D->rmap->n != D->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of D %D do not equal local columns %D",D->rmap->n,D->cmap->n);
-    if (C->rmap->n != D->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of C %D do not equal local rows D %D",C->rmap->n,D->rmap->n);
+  if (!S->assembled) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Use MatSchurComplementSetSubMatrices() for a new matrix");
+  PetscValidHeaderSpecific(A00,MAT_CLASSID,1);
+  PetscValidHeaderSpecific(A01,MAT_CLASSID,2);
+  PetscValidHeaderSpecific(A10,MAT_CLASSID,3);
+  PetscCheckSameComm(A00,1,Ap00,2);
+  PetscCheckSameComm(A00,1,A01,3);
+  PetscCheckSameComm(A00,1,A10,4);
+  if (A00->rmap->n != A00->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of A00 %D do not equal local columns %D",A00->rmap->n,A00->cmap->n);
+  if (A00->rmap->n != Ap00->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of A00 %D do not equal local rows of Ap00 %D",A00->rmap->n,Ap00->rmap->n);
+  if (Ap00->rmap->n != Ap00->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of Ap00 %D do not equal local columns %D",Ap00->rmap->n,Ap00->cmap->n);
+  if (A00->cmap->n != A01->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local columns of A00 %D do not equal local rows of A01 %D",A00->cmap->n,A01->rmap->n);
+  if (A10->cmap->n != A00->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local columns of A10 %D do not equal local rows of A00 %D",A10->cmap->n,A00->rmap->n);
+  if (A11) {
+    PetscValidHeaderSpecific(A11,MAT_CLASSID,5);
+    PetscCheckSameComm(A00,1,A11,5);
+    if (A11->rmap->n != A11->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of A11 %D do not equal local columns %D",A11->rmap->n,A11->cmap->n);
+    if (A10->rmap->n != A11->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of A10 %D do not equal local rows A11 %D",A10->rmap->n,A11->rmap->n);
   }
 
-  ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr);
-  ierr = PetscObjectReference((PetscObject)Ap);CHKERRQ(ierr);
-  ierr = PetscObjectReference((PetscObject)B);CHKERRQ(ierr);
-  ierr = PetscObjectReference((PetscObject)C);CHKERRQ(ierr);
-  if (D) {
-    ierr = PetscObjectReference((PetscObject)D);CHKERRQ(ierr);
+  ierr = PetscObjectReference((PetscObject)A00);CHKERRQ(ierr);
+  ierr = PetscObjectReference((PetscObject)Ap00);CHKERRQ(ierr);
+  ierr = PetscObjectReference((PetscObject)A01);CHKERRQ(ierr);
+  ierr = PetscObjectReference((PetscObject)A10);CHKERRQ(ierr);
+  if (A11) {
+    ierr = PetscObjectReference((PetscObject)A11);CHKERRQ(ierr);
   }
 
   ierr = MatDestroy(&Na->A);CHKERRQ(ierr);
   ierr = MatDestroy(&Na->C);CHKERRQ(ierr);
   ierr = MatDestroy(&Na->D);CHKERRQ(ierr);
 
-  Na->A  = A;
-  Na->Ap = Ap;
-  Na->B  = B;
-  Na->C  = C;
-  Na->D  = D;
+  Na->A  = A00;
+  Na->Ap = Ap00;
+  Na->B  = A01;
+  Na->C  = A10;
+  Na->D  = A11;
 
-  ierr = KSPSetOperators(Na->ksp,A,Ap,str);CHKERRQ(ierr);
+  ierr = KSPSetOperators(Na->ksp,A00,Ap00,str);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
 
 #undef __FUNCT__
-#define __FUNCT__ "MatSchurComplementGetSubmatrices"
+#define __FUNCT__ "MatSchurComplementGetSubMatrices"
 /*@C
-  MatSchurComplementGetSubmatrices - Get the individual submatrices in the Schur complement
+  MatSchurComplementGetSubMatrices - Get the individual submatrices in the Schur complement
 
   Collective on Mat
 
-  Input Parameters:
-+ N - the matrix obtained with MatCreateSchurComplement()
-- A,B,C,D  - the four parts of the original matrix (D is optional)
+  Input Parameter:
+. S                - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
 
-  Note:
-  D is optional, and thus can be NULL
+  Output Paramters:
++ A00,A01,A10,A11  - the four parts of the original matrix A = [A00 A01; A10 A11] (A11 is optional)
+- Ap00             - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A^{-1}.
+
+  Note: A11 is optional, and thus can be NULL.  The submatrices are not increfed before they are returned and should not be modified or destroyed.
 
   Level: intermediate
 
-.seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatCreateSchurComplement(), MatSchurComplementUpdate()
+.seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatCreateSchurComplement(), MatSchurComplementUpdateSubMatrices()
 @*/
-PetscErrorCode  MatSchurComplementGetSubmatrices(Mat N,Mat *A,Mat *Ap,Mat *B,Mat *C,Mat *D)
+PetscErrorCode  MatSchurComplementGetSubMatrices(Mat S,Mat *A00,Mat *Ap00,Mat *A01,Mat *A10,Mat *A11)
 {
-  Mat_SchurComplement *Na = (Mat_SchurComplement*) N->data;
+  Mat_SchurComplement *Na = (Mat_SchurComplement*) S->data;
   PetscErrorCode      ierr;
   PetscBool           flg;
 
   PetscFunctionBegin;
-  PetscValidHeaderSpecific(N,MAT_CLASSID,1);
-  ierr = PetscObjectTypeCompare((PetscObject)N,MATSCHURCOMPLEMENT,&flg);CHKERRQ(ierr);
+  PetscValidHeaderSpecific(S,MAT_CLASSID,1);
+  ierr = PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&flg);CHKERRQ(ierr);
   if (flg) {
-    if (A) *A = Na->A;
-    if (Ap) *Ap = Na->Ap;
-    if (B) *B = Na->B;
-    if (C) *C = Na->C;
-    if (D) *D = Na->D;
+    if (A00) *A00 = Na->A;
+    if (Ap00) *Ap00 = Na->Ap;
+    if (A01) *A01 = Na->B;
+    if (A10) *A10 = Na->C;
+    if (A11) *A11 = Na->D;
   } else {
-    if (A) *A = 0;
-    if (Ap) *Ap = 0;
-    if (B) *B = 0;
-    if (C) *C = 0;
-    if (D) *D = 0;
+    if (A00) *A00 = 0;
+    if (Ap00) *Ap00 = 0;
+    if (A01) *A01 = 0;
+    if (A10) *A10 = 0;
+    if (A11) *A11 = 0;
   }
   PetscFunctionReturn(0);
 }
 #undef __FUNCT__
 #define __FUNCT__ "MatGetSchurComplement_Basic"
 /* Developer Notes: This should be implemented with a MatCreate_SchurComplement() as that is the standard design for new Mat classes. */
-PetscErrorCode MatGetSchurComplement_Basic(Mat mat,IS isrow0,IS iscol0,IS isrow1,IS iscol1,MatReuse mreuse,Mat *newmat,MatReuse preuse,Mat *newpmat)
+PetscErrorCode MatGetSchurComplement_Basic(Mat mat,IS isrow0,IS iscol0,IS isrow1,IS iscol1,MatReuse mreuse,Mat *newmat,MatSchurComplementAinvType ainvtype, MatReuse preuse,Mat *newpmat)
 {
   PetscErrorCode ierr;
   Mat            A=0,Ap=0,B=0,C=0,D=0;
   if (mreuse != MAT_IGNORE_MATRIX) {
     /* Use MatSchurComplement */
     if (mreuse == MAT_REUSE_MATRIX) {
-      ierr = MatSchurComplementGetSubmatrices(*newmat,&A,&Ap,&B,&C,&D);CHKERRQ(ierr);
+      ierr = MatSchurComplementGetSubMatrices(*newmat,&A,&Ap,&B,&C,&D);CHKERRQ(ierr);
       if (!A || !Ap || !B || !C) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Attempting to reuse matrix but Schur complement matrices unset");
       if (A != Ap) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Preconditioning matrix does not match operator");
       ierr = MatDestroy(&Ap);CHKERRQ(ierr); /* get rid of extra reference */
       ierr = MatCreateSchurComplement(A,A,B,C,D,newmat);CHKERRQ(ierr);
       break;
     case MAT_REUSE_MATRIX:
-      ierr = MatSchurComplementUpdate(*newmat,A,A,B,C,D,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
+      ierr = MatSchurComplementUpdateSubMatrices(*newmat,A,A,B,C,D,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
       break;
     default:
       SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Unrecognized value of mreuse");
     }
   }
   if (preuse != MAT_IGNORE_MATRIX) {
-    /* Use the diagonal part of A to form D - C inv(diag(A)) B */
-    Mat         Ad,AdB,S;
-    Vec         diag;
-    PetscInt    i,m,n,mstart,mend;
-    PetscScalar *x;
-
-    /* We could compose these with newpmat so that the matrices can be reused. */
-    if (!A) {ierr = MatGetSubMatrix(mat,isrow0,iscol0,MAT_INITIAL_MATRIX,&A);CHKERRQ(ierr);}
-    if (!B) {ierr = MatGetSubMatrix(mat,isrow0,iscol1,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr);}
-    if (!C) {ierr = MatGetSubMatrix(mat,isrow1,iscol0,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr);}
-    if (!D) {ierr = MatGetSubMatrix(mat,isrow1,iscol1,MAT_INITIAL_MATRIX,&D);CHKERRQ(ierr);}
-
-    ierr = MatGetVecs(A,&diag,NULL);CHKERRQ(ierr);
-    ierr = MatGetDiagonal(A,diag);CHKERRQ(ierr);
-    ierr = VecReciprocal(diag);CHKERRQ(ierr);
-    ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr);
-    /* We need to compute S = D - C inv(diag(A)) B.  For row-oriented formats, it is easy to scale the rows of B and
-     * for column-oriented formats the columns of C can be scaled.  Would skip creating a silly diagonal matrix. */
-    ierr = MatCreate(PetscObjectComm((PetscObject)A),&Ad);CHKERRQ(ierr);
-    ierr = MatSetSizes(Ad,m,n,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
-    ierr = MatSetOptionsPrefix(Ad,((PetscObject)mat)->prefix);CHKERRQ(ierr);
-    ierr = MatAppendOptionsPrefix(Ad,"diag_");CHKERRQ(ierr);
-    ierr = MatSetFromOptions(Ad);CHKERRQ(ierr);
-    ierr = MatSeqAIJSetPreallocation(Ad,1,NULL);CHKERRQ(ierr);
-    ierr = MatMPIAIJSetPreallocation(Ad,1,NULL,0,NULL);CHKERRQ(ierr);
-    ierr = MatGetOwnershipRange(Ad,&mstart,&mend);CHKERRQ(ierr);
-    ierr = VecGetArray(diag,&x);CHKERRQ(ierr);
-    for (i=mstart; i<mend; i++) {
-      ierr = MatSetValue(Ad,i,i,x[i-mstart],INSERT_VALUES);CHKERRQ(ierr);
-    }
-    ierr = VecRestoreArray(diag,&x);CHKERRQ(ierr);
-    ierr = MatAssemblyBegin(Ad,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
-    ierr = MatAssemblyEnd(Ad,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
-    ierr = VecDestroy(&diag);CHKERRQ(ierr);
-
-    ierr     = MatMatMult(Ad,B,MAT_INITIAL_MATRIX,1,&AdB);CHKERRQ(ierr);
-    S        = (preuse == MAT_REUSE_MATRIX) ? *newpmat : (Mat)0;
-    ierr     = MatMatMult(C,AdB,preuse,PETSC_DEFAULT,&S);CHKERRQ(ierr);
-    ierr     = MatAYPX(S,-1,D,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
-    *newpmat = S;
-    ierr     = MatDestroy(&Ad);CHKERRQ(ierr);
-    ierr     = MatDestroy(&AdB);CHKERRQ(ierr);
+    ierr = MatCreateSchurComplementPmat(A,B,C,D,ainvtype,preuse,newpmat);CHKERRQ(ierr);
   }
   ierr = MatDestroy(&A);CHKERRQ(ierr);
   ierr = MatDestroy(&B);CHKERRQ(ierr);
     Collective on Mat
 
     Input Parameters:
-+   mat - Matrix in which the complement is to be taken
++   A      - matrix in which the complement is to be taken
 .   isrow0 - rows to eliminate
 .   iscol0 - columns to eliminate, (isrow0,iscol0) should be square and nonsingular
 .   isrow1 - rows in which the Schur complement is formed
 .   iscol1 - columns in which the Schur complement is formed
 .   mreuse - MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX, use MAT_IGNORE_MATRIX to put nothing in newmat
+.   plump  - the type of approximation used for the inverse of the (0,0) block used in forming Sp:
+$                        MAT_SCHUR_COMPLEMENT_AINV_DIAG | MAT_SCHUR_COMPLEMENT_AINV_LUMP
 -   preuse - MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX, use MAT_IGNORE_MATRIX to put nothing in newpmat
 
     Output Parameters:
-+   newmat - exact Schur complement, often of type MATSCHURCOMPLEMENT which is difficult to use for preconditioning
--   newpmat - approximate Schur complement suitable for preconditioning
++   S      - exact Schur complement, often of type MATSCHURCOMPLEMENT which is difficult to use for preconditioning
+-   Sp     - approximate Schur complement suitable for preconditioning
 
     Note:
     Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
-    application-specific information.  The default for assembled matrices is to use the diagonal of the (0,0) block
-    which will rarely produce a scalable algorithm.
+    application-specific information.  The default for assembled matrices is to use the inverse of the diagonal of
+    the (0,0) block A00 in place of A00^{-1}. This rarely produce a scalable algorithm. Optionally, A00 can be lumped
+    before forming inv(diag(A00)).
 
     Sometimes users would like to provide problem-specific data in the Schur complement, usually only for special row
     and column index sets.  In that case, the user should call PetscObjectComposeFunction() to set
 
     Concepts: matrices^submatrices
 
-.seealso: MatGetSubMatrix(), PCFIELDSPLIT, MatCreateSchurComplement()
+.seealso: MatGetSubMatrix(), PCFIELDSPLIT, MatCreateSchurComplement(), MatSchurComplementAinvType
 @*/
-PetscErrorCode  MatGetSchurComplement(Mat mat,IS isrow0,IS iscol0,IS isrow1,IS iscol1,MatReuse mreuse,Mat *newmat,MatReuse preuse,Mat *newpmat)
+PetscErrorCode  MatGetSchurComplement(Mat A,IS isrow0,IS iscol0,IS isrow1,IS iscol1,MatReuse mreuse,Mat *S,MatSchurComplementAinvType ainvtype,MatReuse preuse,Mat *Sp)
 {
   PetscErrorCode ierr,(*f)(Mat,IS,IS,IS,IS,MatReuse,Mat*,MatReuse,Mat*);
 
   PetscFunctionBegin;
-  PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
+  PetscValidHeaderSpecific(A,MAT_CLASSID,1);
   PetscValidHeaderSpecific(isrow0,IS_CLASSID,2);
   PetscValidHeaderSpecific(iscol0,IS_CLASSID,3);
   PetscValidHeaderSpecific(isrow1,IS_CLASSID,4);
   PetscValidHeaderSpecific(iscol1,IS_CLASSID,5);
-  if (mreuse == MAT_REUSE_MATRIX) PetscValidHeaderSpecific(*newmat,MAT_CLASSID,7);
-  if (preuse == MAT_REUSE_MATRIX) PetscValidHeaderSpecific(*newpmat,MAT_CLASSID,9);
-  PetscValidType(mat,1);
-  if (mat->factortype) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
+  if (mreuse == MAT_REUSE_MATRIX) PetscValidHeaderSpecific(*S,MAT_CLASSID,7);
+  if (preuse == MAT_REUSE_MATRIX) PetscValidHeaderSpecific(*Sp,MAT_CLASSID,9);
+  PetscValidType(A,1);
+  if (A->factortype) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
+
+  ierr = PetscObjectQueryFunction((PetscObject)S,"MatGetSchurComplement_C",&f);CHKERRQ(ierr);
+  if (f) {
+    ierr = (*f)(A,isrow0,iscol0,isrow1,iscol1,mreuse,S,preuse,Sp);CHKERRQ(ierr);
+  } else {
+    ierr = MatGetSchurComplement_Basic(A,isrow0,iscol0,isrow1,iscol1,mreuse,S,ainvtype,preuse,Sp);CHKERRQ(ierr);
+  }
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "MatSchurComplementSetAinvType"
+/*@
+    MatSchurComplementSetAinvType - set the type of approximation used for the inverse of the (0,0) block used in forming Sp in MatSchurComplementGetPmat()
+
+    Not collective.
+
+    Input Parameters:
++   S        - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
+-   ainvtype - type of approximation used to form A00inv from A00 when assembling Sp = A11 - A10 A00inv A01:
+$                      MAT_SCHUR_COMPLEMENT_AINV_DIAG or MAT_SCHUR_COMPLEMENT_AINV_LUMP
+
+    Options database:
+    -mat_schur_complement_ainv_type diag | lump
+
+    Note:
+    Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
+    application-specific information.  The default for assembled matrices is to use the inverse of the diagonal of
+    the (0,0) block A00 in place of A00^{-1}. This rarely produce a scalable algorithm. Optionally, A00 can be lumped
+    before forming inv(diag(A00)).
+
+    Level: advanced
+
+    Concepts: matrices^submatrices
+
+.seealso: MatSchurComplementAinvType, MatCreateSchurComplement(), MatGetSchurComplement(), MatSchurComplementGetPmat(), MatSchurComplementGetAinvType()
+@*/
+PetscErrorCode  MatSchurComplementSetAinvType(Mat S,MatSchurComplementAinvType ainvtype)
+{
+  PetscErrorCode      ierr;
+  const char*         t;
+  PetscBool           isschur;
+  Mat_SchurComplement *schur;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(S,MAT_CLASSID,1);
+  ierr = PetscObjectGetType((PetscObject)S,&t);CHKERRQ(ierr);
+  ierr = PetscStrcmp(t,MATSCHURCOMPLEMENT,&isschur);CHKERRQ(ierr);
+  if (!isschur) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected Mat of type MATSCHURCOMPLEMENT, got %s instead",t);
+  schur = (Mat_SchurComplement*)S->data;
+  if (ainvtype != MAT_SCHUR_COMPLEMENT_AINV_DIAG && ainvtype != MAT_SCHUR_COMPLEMENT_AINV_LUMP) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown MatSchurComplementAinvType: %D",ainvtype);
+  schur->ainvtype = ainvtype;
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "MatSchurComplementGetAinvType"
+/*@
+    MatSchurComplementGetAinvType - get the type of approximation for the inverse of the (0,0) block used in forming Sp in MatSchurComplementGetPmat()
+
+    Not collective.
+
+    Input Parameter:
+.   S      - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
+
+    Output Parameter:
+.   ainvtype - type of approximation used to form A00inv from A00 when assembling Sp = A11 - A10 A00inv A01:
+$                      MAT_SCHUR_COMPLEMENT_AINV_DIAG or MAT_SCHUR_COMPLEMENT_AINV_LUMP
+
+    Note:
+    Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
+    application-specific information.  The default for assembled matrices is to use the inverse of the diagonal of
+    the (0,0) block A00 in place of A00^{-1}. This rarely produce a scalable algorithm. Optionally, A00 can be lumped
+    before forming inv(diag(A00)).
+
+    Level: advanced
+
+    Concepts: matrices^submatrices
+
+.seealso: MatSchurComplementAinvType, MatCreateSchurComplement(), MatGetSchurComplement(), MatSchurComplementGetPmat(), MatSchurComplementSetAinvType()
+@*/
+PetscErrorCode  MatSchurComplementGetAinvType(Mat S,MatSchurComplementAinvType *ainvtype)
+{
+  PetscErrorCode      ierr;
+  const char*         t;
+  PetscBool           isschur;
+  Mat_SchurComplement *schur;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(S,MAT_CLASSID,1);
+  ierr = PetscObjectGetType((PetscObject)S,&t);CHKERRQ(ierr);
+  ierr = PetscStrcmp(t,MATSCHURCOMPLEMENT,&isschur);CHKERRQ(ierr);
+  if (!isschur) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected Mat of type MATSCHURCOMPLEMENT, got %s instead",t);
+  schur = (Mat_SchurComplement*)S->data;
+  if (ainvtype) *ainvtype = schur->ainvtype;
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "MatCreateSchurComplementPmat"
+/*@
+    MatCreateSchurComplementPmat - create a preconditioning matrix for the Schur complement by assembling Sp = A11 - A10 inv(diag(A00)) A01
+
+    Collective on Mat
+
+    Input Parameters:
++   A00,A01,A10,A11      - the four parts of the original matrix A = [A00 A01; A10 A11] (A01,A10, and A11 are optional, implying zero matrices)
+.   ainvtype             - type of approximation for inv(A00) used when forming Sp = A11 - A10 inv(A00) A01
+-   preuse               - MAT_INITIAL_MATRIX for a new Sp, or MAT_REUSE_MATRIX to reuse an existing Sp, or MAT_IGNORE_MATRIX to put nothing in Sp
+
+    Output Parameter:
+-   Spmat                - approximate Schur complement suitable for preconditioning S = A11 - A10 inv(diag(A00)) A01
+
+    Note:
+    Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
+    application-specific information.  The default for assembled matrices is to use the inverse of the diagonal of
+    the (0,0) block A00 in place of A00^{-1}. This rarely produce a scalable algorithm. Optionally, A00 can be lumped
+    before forming inv(diag(A00)).
+
+    Level: advanced
+
+    Concepts: matrices^submatrices
+
+.seealso: MatCreateSchurComplement(), MatGetSchurComplement(), MatSchurComplementGetPmat(), MatSchurComplementAinvType
+@*/
+PetscErrorCode  MatCreateSchurComplementPmat(Mat A00,Mat A01,Mat A10,Mat A11,MatSchurComplementAinvType ainvtype,MatReuse preuse,Mat *Spmat)
+{
+
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  /* Use an appropriate approximate inverse of A00 to form A11 - A10 inv(diag(A00)) A01; a NULL A01, A10 or A11 indicates a zero matrix. */
+  /* TODO: Perhaps should create an appropriately-sized zero matrix of the same type as A00? */
+  if ((!A01 || !A10) & !A11) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Cannot assemble Spmat: A01, A10 and A11 are all NULL.");
+
+  if (!A01 || !A10) {
+    if (!preuse) {
+      ierr = MatDuplicate(A11,MAT_COPY_VALUES,Spmat);CHKERRQ(ierr);
+    } else {
+      ierr = MatCopy(A11,*Spmat,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
+    }
+
+  } else {
+    Mat         AdB,Sp;
+    Vec         diag;
+
+    ierr = MatGetVecs(A00,&diag,NULL);CHKERRQ(ierr);
+    if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_DIAG) {
+      ierr = MatGetRowSum(A00,diag);CHKERRQ(ierr);
+    } else if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_LUMP) {
+      ierr = MatGetDiagonal(A00,diag);CHKERRQ(ierr);
+    } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown MatSchurComplementAinvType: %D", ainvtype);
+    ierr = VecReciprocal(diag);CHKERRQ(ierr);
+    ierr = MatDuplicate(A01,MAT_COPY_VALUES,&AdB);CHKERRQ(ierr);
+    ierr = MatDiagonalScale(AdB,diag,NULL);CHKERRQ(ierr);
+    ierr = VecDestroy(&diag);CHKERRQ(ierr);
+    Sp       = (preuse == MAT_REUSE_MATRIX) ? *Spmat : (Mat)0;
+    ierr     = MatMatMult(A10,AdB,preuse,PETSC_DEFAULT,&Sp);CHKERRQ(ierr);
+    if (!A11) {
+      ierr = MatScale(Sp,-1.0);CHKERRQ(ierr);
+    } else {
+      ierr     = MatAYPX(Sp,-1,A11,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
+    }
+    *Spmat    = Sp;
+    ierr     = MatDestroy(&AdB);CHKERRQ(ierr);
+  }
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "MatSchurComplementGetPmat_Basic"
+PetscErrorCode  MatSchurComplementGetPmat_Basic(Mat S,MatReuse preuse,Mat *Spmat)
+{
+  Mat A,B,C,D;
+  Mat_SchurComplement *schur = (Mat_SchurComplement *)S->data;
+  PetscErrorCode      ierr;
+
+  PetscFunctionBegin;
+  if (preuse == MAT_IGNORE_MATRIX) PetscFunctionReturn(0);
+
+  ierr = MatSchurComplementGetSubMatrices(S,&A,NULL,&B,&C,&D);CHKERRQ(ierr);
+  if (!A) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Schur complement component matrices unset");
+  ierr = MatCreateSchurComplementPmat(A,B,C,D,schur->ainvtype,preuse,Spmat);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "MatSchurComplementGetPmat"
+/*@
+    MatSchurComplementGetPmat - Obtain a preconditioning matrix for the Schur complement by assembling Sp = A11 - A10 inv(diag(A00)) A01
+
+    Collective on Mat
+
+    Input Parameters:
++   S      - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
+-   preuse - MAT_INITIAL_MATRIX for a new Sp, or MAT_REUSE_MATRIX to reuse an existing Sp, or MAT_IGNORE_MATRIX to put nothing in Sp
+
+    Output Parameter:
+-   Sp     - approximate Schur complement suitable for preconditioning S = A11 - A10 inv(diag(A00)) A01
+
+    Note:
+    Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
+    application-specific information.  The default for assembled matrices is to use the inverse of the diagonal of
+    the (0,0) block A00 in place of A00^{-1}. This rarely produce a scalable algorithm. Optionally, A00 can be lumped
+    before forming inv(diag(A00)).
+
+    Sometimes users would like to provide problem-specific data in the Schur complement, usually only
+    for special row and column index sets.  In that case, the user should call PetscObjectComposeFunction() to set
+    "MatSchurComplementGetPmat_C" to their function.  If their function needs to fall back to the default implementation,
+    it should call MatSchurComplementGetPmat_Basic().
+
+    Level: advanced
+
+    Concepts: matrices^submatrices
+
+.seealso: MatGetSubMatrix(), PCFIELDSPLIT, MatGetSchurComplement(), MatCreateSchurComplement(), MatSchurComplementSetAinvType()
+@*/
+PetscErrorCode  MatSchurComplementGetPmat(Mat S,MatReuse preuse,Mat *Sp)
+{
+  PetscErrorCode ierr,(*f)(Mat,MatReuse,Mat*);
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(S,MAT_CLASSID,1);
+  if (preuse == MAT_REUSE_MATRIX) PetscValidHeaderSpecific(*Sp,MAT_CLASSID,3);
+  PetscValidType(S,1);
+  if (S->factortype) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
 
-  ierr = PetscObjectQueryFunction((PetscObject)mat,"MatGetSchurComplement_C",&f);CHKERRQ(ierr);
+  ierr = PetscObjectQueryFunction((PetscObject)S,"MatSchurComplementGetPmat_C",&f);CHKERRQ(ierr);
   if (f) {
-    ierr = (*f)(mat,isrow0,iscol0,isrow1,iscol1,mreuse,newmat,preuse,newpmat);CHKERRQ(ierr);
+    ierr = (*f)(S,preuse,Sp);CHKERRQ(ierr);
   } else {
-    ierr = MatGetSchurComplement_Basic(mat,isrow0,iscol0,isrow1,iscol1,mreuse,newmat,preuse,newpmat);CHKERRQ(ierr);
+    ierr = MatSchurComplementGetPmat_Basic(S,preuse,Sp);CHKERRQ(ierr);
   }
   PetscFunctionReturn(0);
 }

File src/ksp/pc/impls/fieldsplit/fieldsplit.c

        http://chess.cs.umd.edu/~elman/papers/tax.pdf
 */
 
-const char *const PCFieldSplitSchurPreTypes[] = {"SELF","A11","USER","FULL","PCFieldSplitSchurPreType","PC_FIELDSPLIT_SCHUR_PRE_",0};
+const char *const PCFieldSplitSchurPreTypes[] = {"SELF","SELFP","A11","USER","FULL","PCFieldSplitSchurPreType","PC_FIELDSPLIT_SCHUR_PRE_",0};
 const char *const PCFieldSplitSchurFactTypes[] = {"DIAG","LOWER","UPPER","FULL","PCFieldSplitSchurFactType","PC_FIELDSPLIT_SCHUR_FACT_",0};
 
 typedef struct _PC_FieldSplitLink *PC_FieldSplitLink;
   Mat                       B;                     /* The (0,1) block */
   Mat                       C;                     /* The (1,0) block */
   Mat                       schur;                 /* The Schur complement S = A11 - A10 A00^{-1} A01, the KSP here, kspinner, is H_1 in [El08] */
+  Mat                       schurp;                /* Assembled approximation to S built by MatSchurComplement to be used as a preconditioning matrix when solving with S */
   Mat                       schur_user;            /* User-provided preconditioning matrix for the Schur complement */
   PCFieldSplitSchurPreType  schurpre;              /* Determines which preconditioning matrix is used for the Schur complement */
   PCFieldSplitSchurFactType schurfactorization;
 {
   switch (jac->schurpre) {
   case PC_FIELDSPLIT_SCHUR_PRE_SELF: return jac->schur;
+  case PC_FIELDSPLIT_SCHUR_PRE_SELFP: return jac->schurp;
   case PC_FIELDSPLIT_SCHUR_PRE_A11: return jac->pmat[1];
   case PC_FIELDSPLIT_SCHUR_PRE_FULL: /* We calculate this and store it in schur_user */
   case PC_FIELDSPLIT_SCHUR_PRE_USER: /* Use a user-provided matrix if it is given, otherwise diagonal block */
     switch (jac->schurpre) {
     case PC_FIELDSPLIT_SCHUR_PRE_SELF:
       ierr = PetscViewerASCIIPrintf(viewer,"  Preconditioner for the Schur complement formed from S itself\n");CHKERRQ(ierr);break;
+    case PC_FIELDSPLIT_SCHUR_PRE_SELFP:
+      ierr = PetscViewerASCIIPrintf(viewer,"  Preconditioner for the Schur complement formed from Sp, an assembled approximation to S, which uses A00's diagonal's inverse\n");CHKERRQ(ierr);break;
     case PC_FIELDSPLIT_SCHUR_PRE_A11:
       ierr = PetscViewerASCIIPrintf(viewer,"  Preconditioner for the Schur complement formed from A11\n");CHKERRQ(ierr);break;
     case PC_FIELDSPLIT_SCHUR_PRE_FULL:
       ierr  = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr);
       ierr  = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);CHKERRQ(ierr);
       ierr  = ISDestroy(&ccis);CHKERRQ(ierr);
-      ierr  = MatSchurComplementUpdate(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1],pc->flag);CHKERRQ(ierr);
+      ierr  = MatSchurComplementUpdateSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1],pc->flag);CHKERRQ(ierr);
+      if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
+	ierr = MatDestroy(&jac->schurp);CHKERRQ(ierr);
+	ierr = MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);CHKERRQ(ierr);
+      }
       if (kspA != kspInner) {
         ierr = KSPSetOperators(kspA,jac->mat[0],jac->pmat[0],pc->flag);CHKERRQ(ierr);
       }
       ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),pc->flag);CHKERRQ(ierr);
     } else {
       const char   *Dprefix;
-      char         schurprefix[256];
+      char         schurprefix[256], schurmatprefix[256];
       char         schurtestoption[256];
       MatNullSpace sp;
       PetscBool    flg;
       /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */
       ierr = MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur);CHKERRQ(ierr);
       ierr = MatSetType(jac->schur,MATSCHURCOMPLEMENT);CHKERRQ(ierr);
-      ierr = MatSchurComplementSet(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);CHKERRQ(ierr);
-
+      ierr = MatSchurComplementSetSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);CHKERRQ(ierr);
+      ierr = PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr);
+      /* Note that the inner KSP is NOT going to inherit this prefix, and if it did, it would be reset just below.  Is that what we want? */
+      ierr = MatSetOptionsPrefix(jac->schur,schurmatprefix);CHKERRQ(ierr);
+      ierr = MatSetFromOptions(jac->schur);CHKERRQ(ierr);
       ierr = MatGetNullSpace(jac->pmat[1], &sp);CHKERRQ(ierr);
       if (sp) {
         ierr = MatSetNullSpace(jac->schur, sp);CHKERRQ(ierr);
         ierr          = PetscObjectReference((PetscObject) jac->head->ksp);CHKERRQ(ierr);
       }
 
+      if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
+	ierr = MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);CHKERRQ(ierr);
+      }
       ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&jac->kspschur);CHKERRQ(ierr);
       ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);CHKERRQ(ierr);
       ierr = PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);CHKERRQ(ierr);
   ierr       = VecDestroy(&jac->w1);CHKERRQ(ierr);
   ierr       = VecDestroy(&jac->w2);CHKERRQ(ierr);
   ierr       = MatDestroy(&jac->schur);CHKERRQ(ierr);
+  ierr       = MatDestroy(&jac->schurp);CHKERRQ(ierr);
   ierr       = MatDestroy(&jac->schur_user);CHKERRQ(ierr);
   ierr       = KSPDestroy(&jac->kspschur);CHKERRQ(ierr);
   ierr       = KSPDestroy(&jac->kspupper);CHKERRQ(ierr);
     Collective on PC
 
     Input Parameters:
-+   pc  - the preconditioner context
-.   ptype - which matrix to use for preconditioning the Schur complement, PC_FIELDSPLIT_SCHUR_PRE_A11 (diag) is default
++   pc      - the preconditioner context
+.   ptype   - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11 (default), PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_PRE_USER
 -   userpre - matrix to use for preconditioning, or NULL
 
     Options Database:
-.     -pc_fieldsplit_schur_precondition <self,user,a11,full> default is a11
+.     -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> default is a11
 
     Notes:
 $    If ptype is
 $        self the preconditioner for the Schur complement is generated from the Schur complement matrix itself:
 $             The only preconditioner that currently works directly with the Schur complement matrix object is the PCLSC
 $             preconditioner
+$        selfp then the preconditioning matrix is an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01
+$             This is only a good preconditioner when diag(A00) is a good preconditioner for A00.
 
      When solving a saddle point problem, where the A11 block is identically zero, using a11 as the ptype only makes sense
     with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and
-    -fieldsplit_1_pc_type lsc which uses the least squares commutator compute a preconditioner for the Schur complement.
+    -fieldsplit_1_pc_type lsc which uses the least squares commutator to compute a preconditioner for the Schur complement.
 
-    Developer Notes: This is a terrible name, gives no good indication of what the function does and should also have Set in
-     the name since it sets a proceedure to use.
+    Developer Notes: This is a terrible name, which gives no good indication of what the function does.
+                     Among other things, it should have 'Set' in the name, since it sets the type of matrix to use.
 
     Level: intermediate
 
 }
 
 #undef __FUNCT__
+#define __FUNCT__ "PCFieldSplitSchurGetS"
+/*@
+    PCFieldSplitSchurGetS -  extract the MatSchurComplement object used by this PC in case it needs to be configured separately
+
+    Not collective
+
+    Input Parameter:
+.   pc      - the preconditioner context
+
+    Output Parameter:
+.   S       - the Schur complement matrix
+
+    Notes:
+    This matrix should not be destroyed using MatDestroy(); rather, use PCFieldSplitSchurRestoreS().
+
+    Level: advanced
+
+.seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSchurPrecondition(), MatSchurComplement, PCFieldSplitSchurRestoreS()
+
+@*/
+PetscErrorCode  PCFieldSplitSchurGetS(PC pc,Mat *S)
+{
+  PetscErrorCode ierr;
+  const char*    t;
+  PetscBool      isfs;
+  PC_FieldSplit  *jac;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(pc,PC_CLASSID,1);
+  ierr = PetscObjectGetType((PetscObject)pc,&t);CHKERRQ(ierr);
+  ierr = PetscStrcmp(t,PCFIELDSPLIT,&isfs);CHKERRQ(ierr);
+  if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
+  jac = (PC_FieldSplit*)pc->data;
+  if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
+  if (S) *S = jac->schur;
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "PCFieldSplitSchurRestoreS"
+/*@
+    PCFieldSplitSchurRestoreS -  restores the MatSchurComplement object used by this PC
+
+    Not collective
+
+    Input Parameters:
++   pc      - the preconditioner context
+.   S       - the Schur complement matrix
+
+    Level: advanced
+
+.seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSchurPrecondition(), MatSchurComplement, PCFieldSplitSchurGetS()
+
+@*/
+PetscErrorCode  PCFieldSplitSchurRestoreS(PC pc,Mat *S)
+{
+  PetscErrorCode ierr;
+  const char*    t;
+  PetscBool      isfs;
+  PC_FieldSplit  *jac;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(pc,PC_CLASSID,1);
+  ierr = PetscObjectGetType((PetscObject)pc,&t);CHKERRQ(ierr);
+  ierr = PetscStrcmp(t,PCFIELDSPLIT,&isfs);CHKERRQ(ierr);
+  if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
+  jac = (PC_FieldSplit*)pc->data;
+  if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
+  if (!S || *S != jac->schur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten");
+  PetscFunctionReturn(0);
+}
+
+
+#undef __FUNCT__
 #define __FUNCT__ "PCFieldSplitSchurPrecondition_FieldSplit"
 static PetscErrorCode  PCFieldSplitSchurPrecondition_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
 {
 
   PetscFunctionBegin;
   jac->schurpre = ptype;
-  if (pre) {
+  if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) {
     ierr            = MatDestroy(&jac->schur_user);CHKERRQ(ierr);
     jac->schur_user = pre;
     ierr            = PetscObjectReference((PetscObject)jac->schur_user);CHKERRQ(ierr);
 -  flg  - boolean indicating whether to use field splits defined by the DM
 
    Options Database Key:
-.  -pc_fieldsplit_dm_splits 
+.  -pc_fieldsplit_dm_splits
 
    Level: Intermediate
 
                               been supplied explicitly by -pc_fieldsplit_%d_fields
 .   -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields)
 .   -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur> - type of relaxation or factorization splitting
-.   -pc_fieldsplit_schur_precondition <self,user,a11,full> - default is a11
+.   -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11
 .   -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero or negative diagonal and uses Schur complement with no preconditioner as the solver
 
 -    Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_
      in providing the SECOND split or 1 if not give). For PCFieldSplitGetKSP() when field number is 0,
      it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default
      A11 is used to construct a preconditioner for S, use PCFieldSplitSchurPrecondition() to turn on or off this
-     option. You can use the preconditioner PCLSC to precondition the Schur complement with -fieldsplit_1_pc_type lsc. The
-     factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above,
+     option. You can use the preconditioner PCLSC to precondition the Schur complement with -fieldsplit_1_pc_type lsc.
+     When option -fieldsplit_schur_precondition selfp is given, an approximation to S is assembled --
+     Sp = A11 - A10 inv(diag(A00)) A01, which has type AIJ and can be used with a variety of preconditioners
+     (e.g., -fieldsplit_1_pc_type asm).
+     The factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above,
      diag gives
 $              ( inv(A00)     0   )
 $              (   0      -ksp(S) )

File src/ksp/pc/impls/lsc/lsc.c

   ierr = KSPSetOptionsPrefix(lsc->kspL,((PetscObject)pc)->prefix);CHKERRQ(ierr);
   ierr = KSPAppendOptionsPrefix(lsc->kspL,"lsc_");CHKERRQ(ierr);
   ierr = KSPSetFromOptions(lsc->kspL);CHKERRQ(ierr);
-  ierr = MatSchurComplementGetSubmatrices(pc->mat,&A,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
+  ierr = MatSchurComplementGetSubMatrices(pc->mat,&A,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
   ierr = MatGetVecs(A,&lsc->x0,&lsc->y0);CHKERRQ(ierr);
   ierr = MatGetVecs(pc->pmat,&lsc->x1,NULL);CHKERRQ(ierr);
   if (lsc->scalediag) {
   ierr = PetscObjectQuery((PetscObject)pc->pmat,"LSC_Lp",(PetscObject*)&Lp);CHKERRQ(ierr);
   if (!Lp) {ierr = PetscObjectQuery((PetscObject)pc->mat,"LSC_Lp",(PetscObject*)&Lp);CHKERRQ(ierr);}
   if (!L) {
-    ierr = MatSchurComplementGetSubmatrices(pc->mat,NULL,NULL,&B,&C,NULL);CHKERRQ(ierr);
+    ierr = MatSchurComplementGetSubMatrices(pc->mat,NULL,NULL,&B,&C,NULL);CHKERRQ(ierr);
     if (!lsc->L) {
       ierr = MatMatMult(C,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&lsc->L);CHKERRQ(ierr);
     } else {
   }
   if (lsc->scale) {
     Mat Ap;
-    ierr = MatSchurComplementGetSubmatrices(pc->mat,NULL,&Ap,NULL,NULL,NULL);CHKERRQ(ierr);
+    ierr = MatSchurComplementGetSubMatrices(pc->mat,NULL,&Ap,NULL,NULL,NULL);CHKERRQ(ierr);
     ierr = MatGetDiagonal(Ap,lsc->scale);CHKERRQ(ierr); /* Should be the mass matrix, but we don't have plumbing for that yet */
     ierr = VecReciprocal(lsc->scale);CHKERRQ(ierr);
   }
   PetscErrorCode ierr;
 
   PetscFunctionBegin;
-  ierr = MatSchurComplementGetSubmatrices(pc->mat,&A,NULL,&B,&C,NULL);CHKERRQ(ierr);
+  ierr = MatSchurComplementGetSubMatrices(pc->mat,&A,NULL,&B,&C,NULL);CHKERRQ(ierr);
   ierr = KSPSolve(lsc->kspL,x,lsc->x1);CHKERRQ(ierr);
   ierr = MatMult(B,lsc->x1,lsc->x0);CHKERRQ(ierr);
   if (lsc->scale) {

File src/snes/impls/multiblock/multiblock.c

       ierr  = ISComplement(ilink->is,rstart,rend,&ccis);CHKERRQ(ierr);
       ierr  = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);CHKERRQ(ierr);
       ierr  = ISDestroy(&ccis);CHKERRQ(ierr);
-      ierr  = MatSchurComplementUpdate(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->pmat[1],pc->flag);CHKERRQ(ierr);
+      ierr  = MatSchurComplementUpdateSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->pmat[1],pc->flag);CHKERRQ(ierr);
       ierr  = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),pc->flag);CHKERRQ(ierr);
 
     } else {