Commits

Stefano Zampini committed 674ae81

PCBDDC: many changes

- Source code split
- Added Reset procedures
- Support for different MatStructure flags

Comments (0)

Files changed (11)

include/petscpc.h

 #if defined(PETSC_HAVE_PCBDDC)
 /* Enum defining how to treat the coarse problem */
 typedef enum {SEQUENTIAL_BDDC,REPLICATED_BDDC,PARALLEL_BDDC,MULTILEVEL_BDDC} CoarseProblemType;
+PETSC_EXTERN PetscErrorCode PCBDDCSetPrimalVerticesLocalIS(PC,IS);
 PETSC_EXTERN PetscErrorCode PCBDDCSetCoarseningRatio(PC,PetscInt);
 PETSC_EXTERN PetscErrorCode PCBDDCSetMaxLevels(PC,PetscInt);
 PETSC_EXTERN PetscErrorCode PCBDDCSetNullSpace(PC,MatNullSpace);

src/ksp/pc/impls/bddc/bddc.c

    change how to deal with the coarse problem (PCBDDCSetCoarseEnvironment):
      - simplify coarse problem structure -> PCBDDC or PCREDUDANT, nothing else -> same comm for all levels?
      - remove coarse enums and allow use of PCBDDCGetCoarseKSP
-     - remove metis dependency -> use MatPartitioning for multilevel -> Assemble serial adjacency in ManageLocalBoundaries?
+     - remove metis dependency -> use MatPartitioning for multilevel -> Assemble serial adjacency in PCBDDCAnalyzeInterface?
    code refactoring:
      - pick up better names for static functions
    change options structure:
    man pages
 */
 
-/* ----------------------------------------------------------------------------------------------------------------------------------------------
+/* ---------------------------------------------------------------------------------------------------------------------------------------------- 
    Implementation of BDDC preconditioner based on:
    C. Dohrmann "An approximate BDDC preconditioner", Numerical Linear Algebra with Applications Volume 14, Issue 2, pages 149-168, March 2007
    ---------------------------------------------------------------------------------------------------------------------------------------------- */
 
-#include "bddc.h" /*I "petscpc.h" I*/
+#include "bddc.h" /*I "petscpc.h" I*/  /* includes for fortran wrappers */
+#include "bddcprivate.h"
 #include <petscblaslapack.h>
+
+/* prototypes for static functions contained in bddc.c */
+static PetscErrorCode PCBDDCSetUseExactDirichlet(PC,PetscBool);
+static PetscErrorCode PCBDDCSetLevel(PC,PetscInt);
+static PetscErrorCode PCBDDCCoarseSetUp(PC);
+static PetscErrorCode PCBDDCSetUpCoarseEnvironment(PC,PetscScalar*);
+
 /* -------------------------------------------------------------------------- */
-#undef __FUNCT__
+#undef __FUNCT__  
 #define __FUNCT__ "PCSetFromOptions_BDDC"
 PetscErrorCode PCSetFromOptions_BDDC(PC pc)
 {
   ierr = PetscOptionsBool("-pc_bddc_constraints_only","Use only constraints in coarse space (i.e. discard vertices)","none",pcbddc->constraints_flag,&pcbddc->constraints_flag,NULL);CHKERRQ(ierr);
   ierr = PetscOptionsBool("-pc_bddc_faces_only"      ,"Use only faces among constraints of coarse space (i.e. discard edges)"         ,"none",pcbddc->faces_flag      ,&pcbddc->faces_flag      ,NULL);CHKERRQ(ierr);
   ierr = PetscOptionsBool("-pc_bddc_edges_only"      ,"Use only edges among constraints of coarse space (i.e. discard faces)"         ,"none",pcbddc->edges_flag      ,&pcbddc->edges_flag      ,NULL);CHKERRQ(ierr);
-
   /* Coarse solver context */
-  static const char * const avail_coarse_problems[] = {"sequential","replicated","parallel","multilevel","CoarseProblemType","PC_BDDC_",0}; /* order of choiches depends on ENUM defined in bddc.h */
+  static const char * const avail_coarse_problems[] = {"sequential","replicated","parallel","multilevel","CoarseProblemType","PC_BDDC_",0}; /*order of choiches depends on ENUM defined in bddc.h */
   ierr = PetscOptionsEnum("-pc_bddc_coarse_problem_type","Set coarse problem type","none",avail_coarse_problems,(PetscEnum)pcbddc->coarse_problem_type,(PetscEnum*)&pcbddc->coarse_problem_type,NULL);CHKERRQ(ierr);
-
   /* Two different application of BDDC to the whole set of dofs, internal and interface */
   ierr = PetscOptionsBool("-pc_bddc_switch_preconditioning_type","Switch between M_2 (default) and M_3 preconditioners (as defined by Dohrmann)","none",pcbddc->inexact_prec_type,&pcbddc->inexact_prec_type,NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsBool("-pc_bddc_use_change_of_basis","Use change of basis approach for primal space","none",pcbddc->usechangeofbasis,&pcbddc->usechangeofbasis,NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsBool("-pc_bddc_use_change_on_faces","Use change of basis approach for face constraints","none",pcbddc->usechangeonfaces,&pcbddc->usechangeonfaces,NULL);CHKERRQ(ierr);
-
-  pcbddc->usechangeonfaces = pcbddc->usechangeonfaces && pcbddc->usechangeofbasis;
-
+  ierr = PetscOptionsBool("-pc_bddc_use_change_of_basis","Use change of basis approach for primal space","none",pcbddc->use_change_of_basis,&pcbddc->use_change_of_basis,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-pc_bddc_use_change_on_faces","Use change of basis approach for face constraints","none",pcbddc->use_change_on_faces,&pcbddc->use_change_on_faces,NULL);CHKERRQ(ierr);
+  if (!pcbddc->use_change_of_basis) {
+    pcbddc->use_change_on_faces = PETSC_FALSE;
+  }
   ierr = PetscOptionsInt("-pc_bddc_coarsening_ratio","Set coarsening ratio used in multilevel coarsening","none",pcbddc->coarsening_ratio,&pcbddc->coarsening_ratio,NULL);CHKERRQ(ierr);
   ierr = PetscOptionsInt("-pc_bddc_max_levels","Set maximum number of levels for multilevel","none",pcbddc->max_levels,&pcbddc->max_levels,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-pc_bddc_use_deluxe_scaling","Use deluxe scaling for BDDC","none",pcbddc->use_deluxe_scaling,&pcbddc->use_deluxe_scaling,NULL);CHKERRQ(ierr);
   ierr = PetscOptionsTail();CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 /* -------------------------------------------------------------------------- */
+#undef __FUNCT__  
+#define __FUNCT__ "PCBDDCSetPrimalVerticesLocalIS_BDDC"
+static PetscErrorCode PCBDDCSetPrimalVerticesLocalIS_BDDC(PC pc, IS PrimalVertices)
+{
+  PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
+  PetscErrorCode ierr;
 
-#undef __FUNCT__
+  PetscFunctionBegin;
+  ierr = ISDestroy(&pcbddc->user_primal_vertices);CHKERRQ(ierr);
+  ierr = PetscObjectReference((PetscObject)PrimalVertices);CHKERRQ(ierr);
+  pcbddc->user_primal_vertices = PrimalVertices; 
+  PetscFunctionReturn(0);
+}
+#undef __FUNCT__  
+#define __FUNCT__ "PCBDDCSetPrimalVerticesLocalIS"
+/*@
+ PCBDDCSetPrimalVerticesLocalIS - Set user defined primal vertices in PCBDDC.
+
+   Not collective
+
+   Input Parameters:
++  pc - the preconditioning context
+-  PrimalVertices - index sets of primal vertices in local numbering
+
+   Level: intermediate
+
+   Notes:
+
+.seealso: PCBDDC
+@*/
+PetscErrorCode PCBDDCSetPrimalVerticesLocalIS(PC pc, IS PrimalVertices)
+{
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(pc,PC_CLASSID,1);
+  PetscValidHeaderSpecific(PrimalVertices,IS_CLASSID,2);
+  ierr = PetscTryMethod(pc,"PCBDDCSetPrimalVerticesLocalIS_C",(PC,IS),(pc,PrimalVertices));CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+/* -------------------------------------------------------------------------- */
+#undef __FUNCT__  
 #define __FUNCT__ "PCBDDCSetCoarseProblemType_BDDC"
 static PetscErrorCode PCBDDCSetCoarseProblemType_BDDC(PC pc, CoarseProblemType CPT)
 {
-  PC_BDDC *pcbddc = (PC_BDDC*)pc->data;
+  PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;
 
   PetscFunctionBegin;
-  pcbddc->coarse_problem_type = CPT;
+  pcbddc->coarse_problem_type = CPT; 
   PetscFunctionReturn(0);
 }
 
-#undef __FUNCT__
+#undef __FUNCT__  
 #define __FUNCT__ "PCBDDCSetCoarseProblemType"
 /*@
  PCBDDCSetCoarseProblemType - Set coarse problem type in PCBDDC.
   PetscFunctionReturn(0);
 }
 /* -------------------------------------------------------------------------- */
-#undef __FUNCT__
+#undef __FUNCT__  
 #define __FUNCT__ "PCBDDCSetCoarseningRatio_BDDC"
 static PetscErrorCode PCBDDCSetCoarseningRatio_BDDC(PC pc,PetscInt k)
 {
-  PC_BDDC *pcbddc = (PC_BDDC*)pc->data;
+  PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;
 
   PetscFunctionBegin;
   pcbddc->coarsening_ratio=k;
   PetscFunctionReturn(0);
 }
 
-#undef __FUNCT__
+#undef __FUNCT__  
 #define __FUNCT__ "PCBDDCSetCoarseningRatio"
 /*@
  PCBDDCSetCoarseningRatio - Set coarsening ratio used in multilevel coarsening
 }
 /* -------------------------------------------------------------------------- */
 
-#undef __FUNCT__
+#undef __FUNCT__  
 #define __FUNCT__ "PCBDDCSetMaxLevels_BDDC"
 static PetscErrorCode PCBDDCSetMaxLevels_BDDC(PC pc,PetscInt max_levels)
 {
-  PC_BDDC *pcbddc = (PC_BDDC*)pc->data;
+  PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;
 
   PetscFunctionBegin;
   pcbddc->max_levels=max_levels;
   PetscFunctionReturn(0);
 }
 
-#undef __FUNCT__
+#undef __FUNCT__  
 #define __FUNCT__ "PCBDDCSetMaxLevels"
 /*@
  PCBDDCSetMaxLevels - Sets the maximum number of levels within the multilevel approach.
 }
 /* -------------------------------------------------------------------------- */
 
-#undef __FUNCT__
+#undef __FUNCT__  
 #define __FUNCT__ "PCBDDCSetNullSpace_BDDC"
 static PetscErrorCode PCBDDCSetNullSpace_BDDC(PC pc,MatNullSpace NullSpace)
 {
-  PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
+  PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;
   PetscErrorCode ierr;
 
   PetscFunctionBegin;
   ierr = PetscObjectReference((PetscObject)NullSpace);CHKERRQ(ierr);
   ierr = MatNullSpaceDestroy(&pcbddc->NullSpace);CHKERRQ(ierr);
-
-  pcbddc->NullSpace = NullSpace;
+  pcbddc->NullSpace=NullSpace;
   PetscFunctionReturn(0);
 }
 
-#undef __FUNCT__
+#undef __FUNCT__  
 #define __FUNCT__ "PCBDDCSetNullSpace"
 /*@
  PCBDDCSetNullSpace - Set NullSpace of global operator of BDDC preconditioned mat.
 
   PetscFunctionBegin;
   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
+  PetscValidHeaderSpecific(NullSpace,MAT_NULLSPACE_CLASSID,2);
   ierr = PetscTryMethod(pc,"PCBDDCSetNullSpace_C",(PC,MatNullSpace),(pc,NullSpace));CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 /* -------------------------------------------------------------------------- */
 
-#undef __FUNCT__
+#undef __FUNCT__  
 #define __FUNCT__ "PCBDDCSetDirichletBoundaries_BDDC"
 static PetscErrorCode PCBDDCSetDirichletBoundaries_BDDC(PC pc,IS DirichletBoundaries)
 {
-  PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
+  PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;
   PetscErrorCode ierr;
 
   PetscFunctionBegin;
   ierr = ISDestroy(&pcbddc->DirichletBoundaries);CHKERRQ(ierr);
   ierr = PetscObjectReference((PetscObject)DirichletBoundaries);CHKERRQ(ierr);
-
-  pcbddc->DirichletBoundaries = DirichletBoundaries;
+  pcbddc->DirichletBoundaries=DirichletBoundaries;
   PetscFunctionReturn(0);
 }
 
-#undef __FUNCT__
+#undef __FUNCT__  
 #define __FUNCT__ "PCBDDCSetDirichletBoundaries"
 /*@
  PCBDDCSetDirichletBoundaries - Set index set defining subdomain part (in local ordering)
 
   PetscFunctionBegin;
   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
+  PetscValidHeaderSpecific(DirichletBoundaries,IS_CLASSID,2);
   ierr = PetscTryMethod(pc,"PCBDDCSetDirichletBoundaries_C",(PC,IS),(pc,DirichletBoundaries));CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 /* -------------------------------------------------------------------------- */
 
-#undef __FUNCT__
+#undef __FUNCT__  
 #define __FUNCT__ "PCBDDCSetNeumannBoundaries_BDDC"
 static PetscErrorCode PCBDDCSetNeumannBoundaries_BDDC(PC pc,IS NeumannBoundaries)
 {
-  PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
+  PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;
   PetscErrorCode ierr;
 
   PetscFunctionBegin;
   ierr = ISDestroy(&pcbddc->NeumannBoundaries);CHKERRQ(ierr);
   ierr = PetscObjectReference((PetscObject)NeumannBoundaries);CHKERRQ(ierr);
-
-  pcbddc->NeumannBoundaries = NeumannBoundaries;
+  pcbddc->NeumannBoundaries=NeumannBoundaries;
   PetscFunctionReturn(0);
 }
 
-#undef __FUNCT__
+#undef __FUNCT__  
 #define __FUNCT__ "PCBDDCSetNeumannBoundaries"
 /*@
  PCBDDCSetNeumannBoundaries - Set index set defining subdomain part (in local ordering)
 
   PetscFunctionBegin;
   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
+  PetscValidHeaderSpecific(NeumannBoundaries,IS_CLASSID,2);
   ierr = PetscTryMethod(pc,"PCBDDCSetNeumannBoundaries_C",(PC,IS),(pc,NeumannBoundaries));CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 /* -------------------------------------------------------------------------- */
 
-#undef __FUNCT__
+#undef __FUNCT__  
 #define __FUNCT__ "PCBDDCGetDirichletBoundaries_BDDC"
 static PetscErrorCode PCBDDCGetDirichletBoundaries_BDDC(PC pc,IS *DirichletBoundaries)
 {
-  PC_BDDC *pcbddc = (PC_BDDC*)pc->data;
+  PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;
 
   PetscFunctionBegin;
   *DirichletBoundaries = pcbddc->DirichletBoundaries;
   PetscFunctionReturn(0);
 }
 
-#undef __FUNCT__
+#undef __FUNCT__  
 #define __FUNCT__ "PCBDDCGetDirichletBoundaries"
 /*@
  PCBDDCGetDirichletBoundaries - Get index set defining subdomain part (in local ordering)
 }
 /* -------------------------------------------------------------------------- */
 
-#undef __FUNCT__
+#undef __FUNCT__  
 #define __FUNCT__ "PCBDDCGetNeumannBoundaries_BDDC"
 static PetscErrorCode PCBDDCGetNeumannBoundaries_BDDC(PC pc,IS *NeumannBoundaries)
 {
-  PC_BDDC *pcbddc = (PC_BDDC*)pc->data;
+  PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;
 
   PetscFunctionBegin;
   *NeumannBoundaries = pcbddc->NeumannBoundaries;
   PetscFunctionReturn(0);
 }
 
-#undef __FUNCT__
+#undef __FUNCT__  
 #define __FUNCT__ "PCBDDCGetNeumannBoundaries"
 /*@
  PCBDDCGetNeumannBoundaries - Get index set defining subdomain part (in local ordering)
 }
 /* -------------------------------------------------------------------------- */
 
-#undef __FUNCT__
+#undef __FUNCT__  
 #define __FUNCT__ "PCBDDCSetLocalAdjacencyGraph_BDDC"
 static PetscErrorCode PCBDDCSetLocalAdjacencyGraph_BDDC(PC pc, PetscInt nvtxs,const PetscInt xadj[],const PetscInt adjncy[], PetscCopyMode copymode)
 {
-  PC_BDDC        *pcbddc  = (PC_BDDC*)pc->data;
-  PCBDDCGraph    mat_graph=pcbddc->mat_graph;
+  PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
+  PCBDDCGraph    mat_graph = pcbddc->mat_graph;
   PetscErrorCode ierr;
 
   PetscFunctionBegin;
-  mat_graph->nvtxs=nvtxs;
-
-  ierr = PetscFree(mat_graph->xadj);CHKERRQ(ierr);
-  ierr = PetscFree(mat_graph->adjncy);CHKERRQ(ierr);
+  /* free old CSR */
+  ierr = PCBDDCGraphResetCSR(mat_graph);CHKERRQ(ierr);
+  /* get CSR into graph structure */
   if (copymode == PETSC_COPY_VALUES) {
-    ierr = PetscMalloc((mat_graph->nvtxs+1)*sizeof(PetscInt),&mat_graph->xadj);CHKERRQ(ierr);
-    ierr = PetscMalloc(xadj[mat_graph->nvtxs]*sizeof(PetscInt),&mat_graph->adjncy);CHKERRQ(ierr);
-    ierr = PetscMemcpy(mat_graph->xadj,xadj,(mat_graph->nvtxs+1)*sizeof(PetscInt));CHKERRQ(ierr);
-    ierr = PetscMemcpy(mat_graph->adjncy,adjncy,xadj[mat_graph->nvtxs]*sizeof(PetscInt));CHKERRQ(ierr);
+    ierr = PetscMalloc((nvtxs+1)*sizeof(PetscInt),&mat_graph->xadj);CHKERRQ(ierr);
+    ierr = PetscMalloc(xadj[nvtxs]*sizeof(PetscInt),&mat_graph->adjncy);CHKERRQ(ierr);
+    ierr = PetscMemcpy(mat_graph->xadj,xadj,(nvtxs+1)*sizeof(PetscInt));CHKERRQ(ierr);
+    ierr = PetscMemcpy(mat_graph->adjncy,adjncy,xadj[nvtxs]*sizeof(PetscInt));CHKERRQ(ierr);
   } else if (copymode == PETSC_OWN_POINTER) {
-    mat_graph->xadj   = (PetscInt*)xadj;
+    mat_graph->xadj = (PetscInt*)xadj;
     mat_graph->adjncy = (PetscInt*)adjncy;
-  } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported copy mode %d\n",copymode);
+  }
   PetscFunctionReturn(0);
 }
 
-#undef __FUNCT__
+#undef __FUNCT__  
 #define __FUNCT__ "PCBDDCSetLocalAdjacencyGraph"
 /*@
  PCBDDCSetLocalAdjacencyGraph - Set CSR graph of local matrix for use of PCBDDC.
 
   PetscFunctionBegin;
   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
-  ierr = MatGetSize(matis->A,&nrows,&ncols);CHKERRQ(ierr);
-  if (nvtxs != nrows) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local adjacency size %d passed in %s differs from local problem size %d!\n",nvtxs,nrows);
-  else {
-    ierr = PetscTryMethod(pc,"PCBDDCSetLocalAdjacencyGraph_C",(PC,PetscInt,const PetscInt[],const PetscInt[],PetscCopyMode),(pc,nvtxs,xadj,adjncy,copymode));CHKERRQ(ierr);
+  PetscValidIntPointer(xadj,3);
+  PetscValidIntPointer(xadj,4);
+  if (copymode != PETSC_COPY_VALUES && copymode != PETSC_OWN_POINTER) {
+    SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported copy mode %d in %s\n",copymode,__FUNCT__);
   }
+  /* pcis info could not be available at this point */
+  ierr = MatGetSize(matis->A,&nrows,&ncols);CHKERRQ(ierr);
+  if (nvtxs != nrows) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local adjacency size %d passed in %s differs from local problem size %d!\n",nvtxs,__FUNCT__,nrows);
+  ierr = PetscTryMethod(pc,"PCBDDCSetLocalAdjacencyGraph_C",(PC,PetscInt,const PetscInt[],const PetscInt[],PetscCopyMode),(pc,nvtxs,xadj,adjncy,copymode));CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 /* -------------------------------------------------------------------------- */
 
-#undef __FUNCT__
+#undef __FUNCT__  
 #define __FUNCT__ "PCBDDCSetDofsSplitting_BDDC"
 static PetscErrorCode PCBDDCSetDofsSplitting_BDDC(PC pc,PetscInt n_is, IS ISForDofs[])
 {
-  PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
-  PetscInt       i;
+  PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;
+  PetscInt i;
   PetscErrorCode ierr;
 
   PetscFunctionBegin;
   /* Destroy ISes if they were already set */
-  for (i=0; i<pcbddc->n_ISForDofs; i++) {
+  for (i=0;i<pcbddc->n_ISForDofs;i++) {
     ierr = ISDestroy(&pcbddc->ISForDofs[i]);CHKERRQ(ierr);
   }
   ierr = PetscFree(pcbddc->ISForDofs);CHKERRQ(ierr);
   /* allocate space then set */
   ierr = PetscMalloc(n_is*sizeof(IS),&pcbddc->ISForDofs);CHKERRQ(ierr);
-  for (i=0; i<n_is; i++) {
+  for (i=0;i<n_is;i++) {
     ierr = PetscObjectReference((PetscObject)ISForDofs[i]);CHKERRQ(ierr);
-
     pcbddc->ISForDofs[i]=ISForDofs[i];
   }
   pcbddc->n_ISForDofs=n_is;
   PetscFunctionReturn(0);
 }
 
-#undef __FUNCT__
+#undef __FUNCT__  
 #define __FUNCT__ "PCBDDCSetDofsSplitting"
 /*@
  PCBDDCSetDofsSplitting - Set index sets defining fields of local mat.
   PetscFunctionReturn(0);
 }
 /* -------------------------------------------------------------------------- */
-#undef __FUNCT__
+#undef __FUNCT__  
 #define __FUNCT__ "PCPreSolve_BDDC"
 /* -------------------------------------------------------------------------- */
 /*
 {
   PetscErrorCode ierr;
   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
-  PC_IS          *pcis   = (PC_IS*)(pc->data);
-  Mat_IS         *matis  = (Mat_IS*)pc->pmat->data;
+  PC_IS          *pcis = (PC_IS*)(pc->data);
+  Mat_IS         *matis = (Mat_IS*)pc->pmat->data;
   Mat            temp_mat;
   IS             dirIS;
   PetscInt       dirsize,i,*is_indices;
 
   PetscFunctionBegin;
   if (x) {
-    ierr     = PetscObjectReference((PetscObject)x);CHKERRQ(ierr);
+    ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr);
     used_vec = x;
   } else {
-    ierr     = PetscObjectReference((PetscObject)pcbddc->temp_solution);CHKERRQ(ierr);
+    ierr = PetscObjectReference((PetscObject)pcbddc->temp_solution);CHKERRQ(ierr);
     used_vec = pcbddc->temp_solution;
-    ierr     = VecSet(used_vec,0.0);CHKERRQ(ierr);
+    ierr = VecSet(used_vec,0.0);CHKERRQ(ierr);
   }
   /* hack into ksp data structure PCPreSolve comes earlier in src/ksp/ksp/interface/itfunc.c */
   if (ksp) {
     ierr = KSPGetInitialGuessNonzero(ksp,&guess_nonzero);CHKERRQ(ierr);
-    if (!guess_nonzero) {
+    if ( !guess_nonzero ) {
       ierr = VecSet(used_vec,0.0);CHKERRQ(ierr);
     }
   }
   /* store the original rhs */
   ierr = VecCopy(rhs,pcbddc->original_rhs);CHKERRQ(ierr);
 
-  /* Take into account zeroed rows -> change rhs and store solution removed */
+  /* Take into account zeroed rows -> change rhs and store solution removed */    
   ierr = MatGetDiagonal(pc->pmat,pcis->vec1_global);CHKERRQ(ierr);
   ierr = VecPointwiseDivide(pcis->vec1_global,rhs,pcis->vec1_global);CHKERRQ(ierr);
   ierr = VecScatterBegin(matis->ctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-  ierr = VecScatterEnd  (matis->ctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
+  ierr = VecScatterEnd(matis->ctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
   ierr = VecScatterBegin(matis->ctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-  ierr = VecScatterEnd  (matis->ctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
+  ierr = VecScatterEnd(matis->ctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
   ierr = PCBDDCGetDirichletBoundaries(pc,&dirIS);CHKERRQ(ierr);
   if (dirIS) {
     ierr = ISGetSize(dirIS,&dirsize);CHKERRQ(ierr);
     ierr = VecGetArray(pcis->vec2_N,&array_diagonal);CHKERRQ(ierr);
     ierr = ISGetIndices(dirIS,(const PetscInt**)&is_indices);CHKERRQ(ierr);
     for (i=0; i<dirsize; i++) array_x[is_indices[i]] = array_diagonal[is_indices[i]];
-
     ierr = ISRestoreIndices(dirIS,(const PetscInt**)&is_indices);CHKERRQ(ierr);
     ierr = VecRestoreArray(pcis->vec2_N,&array_diagonal);CHKERRQ(ierr);
     ierr = VecRestoreArray(pcis->vec1_N,&array_x);CHKERRQ(ierr);
   }
   ierr = VecScatterBegin(matis->ctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
-  ierr = VecScatterEnd  (matis->ctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
+  ierr = VecScatterEnd(matis->ctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
 
   /* remove the computed solution from the rhs */
   ierr = VecScale(used_vec,-1.0);CHKERRQ(ierr);
   }
 
   /* rhs change of basis */
-  if (pcbddc->usechangeofbasis) {
+  if (pcbddc->use_change_of_basis) {
     /* swap pointers for local matrices */
-    temp_mat          = matis->A;
-    matis->A          = pcbddc->local_mat;
+    temp_mat = matis->A;
+    matis->A = pcbddc->local_mat;
     pcbddc->local_mat = temp_mat;
     /* Get local rhs and apply transformation of basis */
     ierr = VecScatterBegin(pcis->global_to_B,rhs,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
     /* put back modified values into the global vec using INSERT_VALUES copy mode */
     ierr = VecScatterBegin(pcis->global_to_B,pcis->vec2_B,rhs,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
     ierr = VecScatterEnd  (pcis->global_to_B,pcis->vec2_B,rhs,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
-    if (ksp && pcbddc->NullSpace) {
-      ierr = MatNullSpaceRemove(pcbddc->NullSpace,used_vec,NULL);CHKERRQ(ierr);
-      ierr = MatNullSpaceRemove(pcbddc->NullSpace,rhs,NULL);CHKERRQ(ierr);
-    }
+  }
+  if (ksp && pcbddc->NullSpace) {
+    ierr = MatNullSpaceRemove(pcbddc->NullSpace,used_vec,NULL);CHKERRQ(ierr);
+    ierr = MatNullSpaceRemove(pcbddc->NullSpace,rhs,NULL);CHKERRQ(ierr);
   }
   ierr = VecDestroy(&used_vec);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 /* -------------------------------------------------------------------------- */
-#undef __FUNCT__
+#undef __FUNCT__  
 #define __FUNCT__ "PCPostSolve_BDDC"
 /* -------------------------------------------------------------------------- */
 /*
   PetscErrorCode ierr;
   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
   PC_IS          *pcis   = (PC_IS*)(pc->data);
-  Mat_IS         *matis  = (Mat_IS*)pc->pmat->data;
+  Mat_IS         *matis = (Mat_IS*)pc->pmat->data;
   Mat            temp_mat;
 
   PetscFunctionBegin;
-  if (pcbddc->usechangeofbasis) {
+  if (pcbddc->use_change_of_basis) {
     /* swap pointers for local matrices */
-    temp_mat          = matis->A;
-    matis->A          = pcbddc->local_mat;
+    temp_mat = matis->A;
+    matis->A = pcbddc->local_mat;
     pcbddc->local_mat = temp_mat;
     /* restore rhs to its original state */
     if (rhs) {
   PetscFunctionReturn(0);
 }
 /* -------------------------------------------------------------------------- */
-#undef __FUNCT__
+#undef __FUNCT__  
 #define __FUNCT__ "PCSetUp_BDDC"
 /* -------------------------------------------------------------------------- */
 /*
    PCSetUp_BDDC - Prepares for the use of the BDDC preconditioner
-                  by setting data structures and options.
+                  by setting data structures and options.   
 
    Input Parameter:
 +  pc - the preconditioner context
 PetscErrorCode PCSetUp_BDDC(PC pc)
 {
   PetscErrorCode ierr;
-  PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
+  PC_BDDC*       pcbddc = (PC_BDDC*)pc->data;
+  MatStructure   flag;
+  PetscBool      computeis,computetopography,computesolvers;
 
   PetscFunctionBegin;
-  if (!pc->setupcalled) {
-    /* For BDDC we need to define a local "Neumann" problem different to that defined in PCISSetup
-       So, we set to pcnone the Neumann problem of pcis in order to avoid unneeded computation
-       Also, we decide to directly build the (same) Dirichlet problem */
-    ierr = PetscOptionsSetValue("-is_localN_pc_type","none");CHKERRQ(ierr);
-    ierr = PetscOptionsSetValue("-is_localD_pc_type","none");CHKERRQ(ierr);
-    /* Set up all the "iterative substructuring" common block */
-
-    ierr = PCISSetUp(pc);CHKERRQ(ierr);
-    /* Get stdout for dbg */
-    if (pcbddc->dbg_flag) {
-      ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc),&pcbddc->dbg_viewer);CHKERRQ(ierr);
-      ierr = PetscViewerASCIISynchronizedAllow(pcbddc->dbg_viewer,PETSC_TRUE);CHKERRQ(ierr);
+  /* the following lines of code should be replaced by a better logic between PCIS, PCNN, PCBDDC and other nonoverlapping preconditioners */
+  /* For BDDC we need to define a local "Neumann" problem different to that defined in PCISSetup
+     So, we set to pcnone the Neumann problem of pcis in order to avoid unneeded computation
+     Also, we decide to directly build the (same) Dirichlet problem */
+  ierr = PetscOptionsSetValue("-is_localN_pc_type","none");CHKERRQ(ierr);
+  ierr = PetscOptionsSetValue("-is_localD_pc_type","none");CHKERRQ(ierr);
+  /* Get stdout for dbg */
+  if (pcbddc->dbg_flag && !pcbddc->dbg_viewer) {
+    ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc),&pcbddc->dbg_viewer);CHKERRQ(ierr);
+    ierr = PetscViewerASCIISynchronizedAllow(pcbddc->dbg_viewer,PETSC_TRUE);CHKERRQ(ierr);
+  }
+  /* first attempt to split work */
+  if (pc->setupcalled) {
+    computeis = PETSC_FALSE;
+    ierr = PCGetOperators(pc,NULL,NULL,&flag);CHKERRQ(ierr);
+    if (flag == SAME_PRECONDITIONER) {
+      computetopography = PETSC_FALSE;
+      computesolvers = PETSC_FALSE;
+    } else if (flag == SAME_NONZERO_PATTERN) {
+      computetopography = PETSC_FALSE;
+      computesolvers = PETSC_TRUE;
+    } else { /* DIFFERENT_NONZERO_PATTERN */
+      computetopography = PETSC_TRUE;
+      computesolvers = PETSC_TRUE;
     }
-    /* Analyze local interface */
-    ierr = PCBDDCManageLocalBoundaries(pc);CHKERRQ(ierr);
-    /* Set up local constraint matrix */
-    ierr = PCBDDCCreateConstraintMatrix(pc);CHKERRQ(ierr);
+  } else {
+    computeis = PETSC_TRUE;
+    computetopography = PETSC_TRUE;
+    computesolvers = PETSC_TRUE;
+  }
+  /* Set up all the "iterative substructuring" common block */
+  if (computeis) {
+    ierr = PCISSetUp(pc);CHKERRQ(ierr);
+  }
+  /* Analyze interface and set up local constraint and change of basis matrices */
+  if (computetopography) {
+    /* reset data */
+    ierr = PCBDDCResetTopography(pc);CHKERRQ(ierr);
+    ierr = PCBDDCAnalyzeInterface(pc);CHKERRQ(ierr);
+    ierr = PCBDDCConstraintsSetUp(pc);CHKERRQ(ierr);
+  }
+  if (computesolvers) {
+    /* reset data */
+    ierr = PCBDDCResetSolvers(pc);CHKERRQ(ierr);
+    ierr = PCBDDCScalingDestroy(pc);CHKERRQ(ierr);
     /* Create coarse and local stuffs used for evaluating action of preconditioner */
     ierr = PCBDDCCoarseSetUp(pc);CHKERRQ(ierr);
+    ierr = PCBDDCScalingSetUp(pc);CHKERRQ(ierr);
   }
   PetscFunctionReturn(0);
 }
 
    Application Interface Routine: PCApply()
  */
-#undef __FUNCT__
+#undef __FUNCT__  
 #define __FUNCT__ "PCApply_BDDC"
 PetscErrorCode PCApply_BDDC(PC pc,Vec r,Vec z)
 {
-  PC_IS             *pcis   = (PC_IS*)(pc->data);
+  PC_IS             *pcis = (PC_IS*)(pc->data);
   PC_BDDC           *pcbddc = (PC_BDDC*)(pc->data);
   PetscErrorCode    ierr;
-  const PetscScalar one   = 1.0;
+  const PetscScalar one = 1.0;
   const PetscScalar m_one = -1.0;
-  const PetscScalar zero  = 0.0;
+  const PetscScalar zero = 0.0;
 
 /* This code is similar to that provided in nn.c for PCNN
    NN interface preconditioner changed to BDDC
     ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr);
     /*
       Assembling right hand side for BDDC operator
-      - vec1_D for the Dirichlet part (if needed, i.e. prec_flag=PETSC_TRUE)
-      - the interface part of the global vector z
+      - pcis->vec1_D for the Dirichlet part (if needed, i.e. prec_flag=PETSC_TRUE)
+      - pcis->vec1_B the interface part of the global vector z
     */
     ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr);
     ierr = MatMult(pcis->A_BI,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr);
     ierr = VecCopy(r,z);CHKERRQ(ierr);
     ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
     ierr = VecScatterEnd  (pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
-    ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-    ierr = VecScatterEnd  (pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
+    ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr);
   } else {
-    ierr = VecScatterBegin(pcis->global_to_B,r,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-    ierr = VecScatterEnd  (pcis->global_to_B,r,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
     ierr = VecSet(pcis->vec1_D,zero);CHKERRQ(ierr);
     ierr = VecSet(pcis->vec2_D,zero);CHKERRQ(ierr);
+    ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr);
   }
 
-  /* Apply partition of unity */
-  ierr = VecPointwiseMult(pcis->vec1_B,pcis->D,pcis->vec1_B);CHKERRQ(ierr);
-
   /* Apply interface preconditioner
      input/output vecs: pcis->vec1_B and pcis->vec1_D */
   ierr = PCBDDCApplyInterfacePreconditioner(pc);CHKERRQ(ierr);
 
-  /* Apply partition of unity and sum boundary values */
-  ierr = VecPointwiseMult(pcis->vec1_B,pcis->D,pcis->vec1_B);CHKERRQ(ierr);
-  ierr = VecSet(z,zero);CHKERRQ(ierr);
-  ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
-  ierr = VecScatterEnd  (pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
+  /* Apply transpose of partition of unity operator */
+  ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr);
 
   /* Second Dirichlet solve and assembling of output */
   ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
   if (pcbddc->inexact_prec_type) { ierr = MatMultAdd(pcis->A_II,pcis->vec1_D,pcis->vec3_D,pcis->vec3_D);CHKERRQ(ierr); }
   ierr = KSPSolve(pcbddc->ksp_D,pcis->vec3_D,pcbddc->vec4_D);CHKERRQ(ierr);
   ierr = VecScale(pcbddc->vec4_D,m_one);CHKERRQ(ierr);
-  if (pcbddc->inexact_prec_type) { ierr = VecAXPY (pcbddc->vec4_D,one,pcis->vec1_D);CHKERRQ(ierr); }
-  ierr = VecAXPY (pcis->vec2_D,one,pcbddc->vec4_D);CHKERRQ(ierr);
+  if (pcbddc->inexact_prec_type) { ierr = VecAXPY (pcbddc->vec4_D,one,pcis->vec1_D);CHKERRQ(ierr); } 
+  ierr = VecAXPY (pcis->vec2_D,one,pcbddc->vec4_D);CHKERRQ(ierr); 
   ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
   ierr = VecScatterEnd  (pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
   PetscFunctionReturn(0);
-
 }
 /* -------------------------------------------------------------------------- */
-#undef __FUNCT__
+
+#undef __FUNCT__  
 #define __FUNCT__ "PCDestroy_BDDC"
 PetscErrorCode PCDestroy_BDDC(PC pc)
 {
   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
-  PetscInt       i;
   PetscErrorCode ierr;
 
   PetscFunctionBegin;
   /* free data created by PCIS */
   ierr = PCISDestroy(pc);CHKERRQ(ierr);
-  /* free BDDC data  */
-  ierr = MatNullSpaceDestroy(&pcbddc->CoarseNullSpace);CHKERRQ(ierr);
-  ierr = MatNullSpaceDestroy(&pcbddc->NullSpace);CHKERRQ(ierr);
-  ierr = VecDestroy(&pcbddc->temp_solution);CHKERRQ(ierr);
-  ierr = VecDestroy(&pcbddc->original_rhs);CHKERRQ(ierr);
-  ierr = MatDestroy(&pcbddc->local_mat);CHKERRQ(ierr);
-  ierr = MatDestroy(&pcbddc->ChangeOfBasisMatrix);CHKERRQ(ierr);
-  ierr = VecDestroy(&pcbddc->coarse_vec);CHKERRQ(ierr);
-  ierr = VecDestroy(&pcbddc->coarse_rhs);CHKERRQ(ierr);
-  ierr = KSPDestroy(&pcbddc->coarse_ksp);CHKERRQ(ierr);
-  ierr = MatDestroy(&pcbddc->coarse_mat);CHKERRQ(ierr);
-  ierr = MatDestroy(&pcbddc->coarse_phi_B);CHKERRQ(ierr);
-  ierr = MatDestroy(&pcbddc->coarse_phi_D);CHKERRQ(ierr);
-  ierr = VecDestroy(&pcbddc->vec1_P);CHKERRQ(ierr);
-  ierr = VecDestroy(&pcbddc->vec1_C);CHKERRQ(ierr);
-  ierr = MatDestroy(&pcbddc->local_auxmat1);CHKERRQ(ierr);
-  ierr = MatDestroy(&pcbddc->local_auxmat2);CHKERRQ(ierr);
-  ierr = VecDestroy(&pcbddc->vec1_R);CHKERRQ(ierr);
-  ierr = VecDestroy(&pcbddc->vec2_R);CHKERRQ(ierr);
-  ierr = VecDestroy(&pcbddc->vec4_D);CHKERRQ(ierr);
-  ierr = VecScatterDestroy(&pcbddc->R_to_B);CHKERRQ(ierr);
-  ierr = VecScatterDestroy(&pcbddc->R_to_D);CHKERRQ(ierr);
-  ierr = VecScatterDestroy(&pcbddc->coarse_loc_to_glob);CHKERRQ(ierr);
-  ierr = KSPDestroy(&pcbddc->ksp_D);CHKERRQ(ierr);
-  ierr = KSPDestroy(&pcbddc->ksp_R);CHKERRQ(ierr);
-  ierr = ISDestroy(&pcbddc->NeumannBoundaries);CHKERRQ(ierr);
-  ierr = ISDestroy(&pcbddc->DirichletBoundaries);CHKERRQ(ierr);
-  ierr = MatDestroy(&pcbddc->ConstraintMatrix);CHKERRQ(ierr);
-  ierr = PetscFree(pcbddc->local_primal_indices);CHKERRQ(ierr);
-  ierr = PetscFree(pcbddc->replicated_local_primal_indices);CHKERRQ(ierr);
-  ierr = PetscFree(pcbddc->replicated_local_primal_values);CHKERRQ(ierr);
-  ierr = PetscFree(pcbddc->local_primal_displacements);CHKERRQ(ierr);
-  ierr = PetscFree(pcbddc->local_primal_sizes);CHKERRQ(ierr);
-  for (i=0; i<pcbddc->n_ISForDofs; i++) {
-    ierr = ISDestroy(&pcbddc->ISForDofs[i]);CHKERRQ(ierr);
-  }
-  ierr = PetscFree(pcbddc->ISForDofs);CHKERRQ(ierr);
-  for (i=0; i<pcbddc->n_ISForFaces; i++) {
-    ierr = ISDestroy(&pcbddc->ISForFaces[i]);CHKERRQ(ierr);
-  }
-  ierr = PetscFree(pcbddc->ISForFaces);CHKERRQ(ierr);
-  for (i=0; i<pcbddc->n_ISForEdges; i++) {
-    ierr = ISDestroy(&pcbddc->ISForEdges[i]);CHKERRQ(ierr);
-  }
-  ierr = PetscFree(pcbddc->ISForEdges);CHKERRQ(ierr);
-  ierr = ISDestroy(&pcbddc->ISForVertices);CHKERRQ(ierr);
-  /* Free graph structure */
-  ierr = PetscFree(pcbddc->mat_graph->xadj);CHKERRQ(ierr);
-  ierr = PetscFree(pcbddc->mat_graph->adjncy);CHKERRQ(ierr);
-  if (pcbddc->mat_graph->nvtxs) {
-    ierr = PetscFree(pcbddc->mat_graph->neighbours_set[0]);CHKERRQ(ierr);
-  }
-  ierr = PetscFree(pcbddc->mat_graph->neighbours_set);CHKERRQ(ierr);
-  ierr = PetscFree4(pcbddc->mat_graph->where,pcbddc->mat_graph->count,pcbddc->mat_graph->cptr,pcbddc->mat_graph->queue);CHKERRQ(ierr);
-  ierr = PetscFree2(pcbddc->mat_graph->which_dof,pcbddc->mat_graph->touched);CHKERRQ(ierr);
-  ierr = PetscFree(pcbddc->mat_graph->where_ncmps);CHKERRQ(ierr);
+  /* free BDDC custom data  */
+  ierr = PCBDDCResetCustomization(pc);CHKERRQ(ierr);
+  /* destroy objects related to topography */
+  ierr = PCBDDCResetTopography(pc);CHKERRQ(ierr);
+  /* free allocated graph structure */
   ierr = PetscFree(pcbddc->mat_graph);CHKERRQ(ierr);
+  /* free data for scaling operator */
+  ierr = PCBDDCScalingDestroy(pc);CHKERRQ(ierr);
+  /* free solvers stuff */
+  ierr = PCBDDCResetSolvers(pc);CHKERRQ(ierr);
   /* remove functions */
+  ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",NULL);CHKERRQ(ierr);
   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",NULL);CHKERRQ(ierr);
   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetMaxLevels_C",NULL);CHKERRQ(ierr);
   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNullSpace_C",NULL);CHKERRQ(ierr);
   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseProblemType_C",NULL);CHKERRQ(ierr);
   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",NULL);CHKERRQ(ierr);
   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",NULL);CHKERRQ(ierr);
-  ierr = PetscObjectComposeFunction((PetscObject)pc,"PCPreSolve_C",NULL);CHKERRQ(ierr);
-  ierr = PetscObjectComposeFunction((PetscObject)pc,"PCPostSolve_C",NULL);CHKERRQ(ierr);
   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",NULL);CHKERRQ(ierr);
   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",NULL);CHKERRQ(ierr);
   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",NULL);CHKERRQ(ierr);
-  /* Free the private data structure that was hanging off the PC */
-  ierr = PetscFree(pcbddc);CHKERRQ(ierr);
+  /* Free the private data structure */
+  ierr = PetscFree(pc->data);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 /* -------------------------------------------------------------------------- */
 
-#undef __FUNCT__
+#undef __FUNCT__  
 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS_BDDC"
 static PetscErrorCode PCBDDCMatFETIDPGetRHS_BDDC(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs)
 {
-  FETIDPMat_ctx  *mat_ctx;
-  PC_IS          * pcis;
-  PC_BDDC        * pcbddc;
+  FETIDPMat_ctx  mat_ctx;
+  PC_IS*         pcis;
+  PC_BDDC*       pcbddc;
   PetscErrorCode ierr;
 
   PetscFunctionBegin;
-  ierr   = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr);
-  pcis   = (PC_IS*)mat_ctx->pc->data;
+  ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr);
+  pcis = (PC_IS*)mat_ctx->pc->data;
   pcbddc = (PC_BDDC*)mat_ctx->pc->data;
 
   /* change of basis for physical rhs if needed
   (*mat_ctx->pc->ops->presolve)(mat_ctx->pc,NULL,standard_rhs,NULL);
   /* store vectors for computation of fetidp final solution */
   ierr = VecScatterBegin(pcis->global_to_D,standard_rhs,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-  ierr = VecScatterEnd  (pcis->global_to_D,standard_rhs,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
+  ierr = VecScatterEnd(pcis->global_to_D,standard_rhs,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
+  /* scale rhs since it should be unassembled : TODO use counter scaling? (also below) */
   ierr = VecScatterBegin(pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-  ierr = VecScatterEnd  (pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-  /* scale rhs since it should be unassembled */
+  ierr = VecScatterEnd(pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
+  /* Apply partition of unity */
   ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr);
+  /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,standard_rhs,mat_ctx->temp_solution_B);CHKERRQ(ierr); */
   if (!pcbddc->inexact_prec_type) {
     /* compute partially subassembled Schur complement right-hand side */
     ierr = KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr);
     ierr = VecAXPY(mat_ctx->temp_solution_B,-1.0,pcis->vec1_B);CHKERRQ(ierr);
     ierr = VecSet(standard_rhs,0.0);CHKERRQ(ierr);
     ierr = VecScatterBegin(pcis->global_to_B,mat_ctx->temp_solution_B,standard_rhs,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
-    ierr = VecScatterEnd  (pcis->global_to_B,mat_ctx->temp_solution_B,standard_rhs,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
+    ierr = VecScatterEnd(pcis->global_to_B,mat_ctx->temp_solution_B,standard_rhs,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
+    /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,standard_rhs,mat_ctx->temp_solution_B);CHKERRQ(ierr); */
     ierr = VecScatterBegin(pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-    ierr = VecScatterEnd  (pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
+    ierr = VecScatterEnd(pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
     ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr);
   }
   /* BDDC rhs */
   PetscFunctionReturn(0);
 }
 
-#undef __FUNCT__
+#undef __FUNCT__  
 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS"
 /*@
  PCBDDCMatFETIDPGetRHS - Get rhs for FETIDP linear system.
 
    Input Parameters:
 +  fetidp_mat   - the FETIDP mat obtained by a call to PCBDDCCreateFETIDPOperators
-+  standard_rhs - the rhs of your linear system
-
++  standard_rhs - the rhs of your linear system 
+  
    Output Parameters:
 +  fetidp_flux_rhs   - the rhs of the FETIDP linear system
 
 @*/
 PetscErrorCode PCBDDCMatFETIDPGetRHS(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs)
 {
-  FETIDPMat_ctx  *mat_ctx;
+  FETIDPMat_ctx  mat_ctx;
   PetscErrorCode ierr;
 
   PetscFunctionBegin;
 }
 /* -------------------------------------------------------------------------- */
 
-#undef __FUNCT__
+#undef __FUNCT__  
 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution_BDDC"
 static PetscErrorCode PCBDDCMatFETIDPGetSolution_BDDC(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol)
 {
-  FETIDPMat_ctx  *mat_ctx;
-  PC_IS          *pcis;
-  PC_BDDC        *pcbddc;
+  FETIDPMat_ctx  mat_ctx;
+  PC_IS*         pcis;
+  PC_BDDC*       pcbddc;
   PetscErrorCode ierr;
 
   PetscFunctionBegin;
-  ierr   = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr);
-  pcis   = (PC_IS*)mat_ctx->pc->data;
+  ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr);
+  pcis = (PC_IS*)mat_ctx->pc->data;
   pcbddc = (PC_BDDC*)mat_ctx->pc->data;
 
   /* apply B_delta^T */
 
 }
 
-#undef __FUNCT__
+#undef __FUNCT__  
 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution"
 /*@
  PCBDDCMatFETIDPGetSolution - Get Solution for FETIDP linear system.
    Input Parameters:
 +  fetidp_mat        - the FETIDP mat obtained by a call to PCBDDCCreateFETIDPOperators
 +  fetidp_flux_sol - the solution of the FETIDP linear system
-
+  
    Output Parameters:
 +  standard_sol      - the solution on the global domain
 
 @*/
 PetscErrorCode PCBDDCMatFETIDPGetSolution(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol)
 {
-  FETIDPMat_ctx  *mat_ctx;
+  FETIDPMat_ctx  mat_ctx;
   PetscErrorCode ierr;
 
   PetscFunctionBegin;
 extern PetscErrorCode PCBDDCDestroyFETIDPMat(Mat);
 extern PetscErrorCode FETIDPPCApply(PC,Vec,Vec);
 extern PetscErrorCode PCBDDCDestroyFETIDPPC(PC);
-#undef __FUNCT__
+
+#undef __FUNCT__  
 #define __FUNCT__ "PCBDDCCreateFETIDPOperators_BDDC"
 static PetscErrorCode PCBDDCCreateFETIDPOperators_BDDC(PC pc, Mat *fetidp_mat, PC *fetidp_pc)
 {
-  FETIDPMat_ctx  *fetidpmat_ctx;
+
+  FETIDPMat_ctx  fetidpmat_ctx;
   Mat            newmat;
-  FETIDPPC_ctx   *fetidppc_ctx;
+  FETIDPPC_ctx   fetidppc_ctx;
   PC             newpc;
   MPI_Comm       comm;
   PetscErrorCode ierr;
   PetscFunctionBegin;
   ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
   /* FETIDP linear matrix */
-  ierr = PCBDDCCreateFETIDPMatContext(pc, &fetidpmat_ctx);CHKERRQ(ierr);
+  ierr = PCBDDCCreateFETIDPMatContext(pc,&fetidpmat_ctx);CHKERRQ(ierr);
   ierr = PCBDDCSetupFETIDPMatContext(fetidpmat_ctx);CHKERRQ(ierr);
   ierr = MatCreateShell(comm,PETSC_DECIDE,PETSC_DECIDE,fetidpmat_ctx->n_lambda,fetidpmat_ctx->n_lambda,fetidpmat_ctx,&newmat);CHKERRQ(ierr);
   ierr = MatShellSetOperation(newmat,MATOP_MULT,(void (*)(void))FETIDPMatMult);CHKERRQ(ierr);
   ierr = MatShellSetOperation(newmat,MATOP_DESTROY,(void (*)(void))PCBDDCDestroyFETIDPMat);CHKERRQ(ierr);
   ierr = MatSetUp(newmat);CHKERRQ(ierr);
   /* FETIDP preconditioner */
-  ierr = PCBDDCCreateFETIDPPCContext(pc, &fetidppc_ctx);CHKERRQ(ierr);
+  ierr = PCBDDCCreateFETIDPPCContext(pc,&fetidppc_ctx);CHKERRQ(ierr);
   ierr = PCBDDCSetupFETIDPPCContext(newmat,fetidppc_ctx);CHKERRQ(ierr);
   ierr = PCCreate(comm,&newpc);CHKERRQ(ierr);
   ierr = PCSetType(newpc,PCSHELL);CHKERRQ(ierr);
   ierr = PCShellSetDestroy(newpc,PCBDDCDestroyFETIDPPC);CHKERRQ(ierr);
   ierr = PCSetOperators(newpc,newmat,newmat,SAME_PRECONDITIONER);CHKERRQ(ierr);
   ierr = PCSetUp(newpc);CHKERRQ(ierr);
-
   /* return pointers for objects created */
-  *fetidp_mat = newmat;
-  *fetidp_pc  = newpc;
+  *fetidp_mat=newmat;
+  *fetidp_pc=newpc;
   PetscFunctionReturn(0);
 }
 
-#undef __FUNCT__
+#undef __FUNCT__  
 #define __FUNCT__ "PCBDDCCreateFETIDPOperators"
 /*@
  PCBDDCCreateFETIDPOperators - Create operators for FETIDP.
 
    Level: intermediate
 
-   Notes: The matrix used with this preconditioner must be of type MATIS
+   Notes: The matrix used with this preconditioner must be of type MATIS 
 
           Unlike more 'conventional' interface preconditioners, this iterates over ALL the
           degrees of freedom, NOT just those on the interface (this allows the use of approximate solvers
 .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,  MATIS
 M*/
 
-#undef __FUNCT__
+#undef __FUNCT__  
 #define __FUNCT__ "PCCreate_BDDC"
 PETSC_EXTERN PetscErrorCode PCCreate_BDDC(PC pc)
 {
-  PetscErrorCode ierr;
-  PC_BDDC        *pcbddc;
-  PCBDDCGraph    mat_graph;
+  PetscErrorCode      ierr;
+  PC_BDDC             *pcbddc;
 
   PetscFunctionBegin;
   /* Creates the private data structure for this preconditioner and attach it to the PC object. */
-  ierr     = PetscNewLog(pc,PC_BDDC,&pcbddc);CHKERRQ(ierr);
-  pc->data = (void*)pcbddc;
+  ierr      = PetscNewLog(pc,PC_BDDC,&pcbddc);CHKERRQ(ierr);
+  pc->data  = (void*)pcbddc;
 
   /* create PCIS data structure */
   ierr = PCISCreate(pc);CHKERRQ(ierr);
 
   /* BDDC specific */
-  pcbddc->CoarseNullSpace                 = 0;
-  pcbddc->NullSpace                       = 0;
-  pcbddc->temp_solution                   = 0;
-  pcbddc->original_rhs                    = 0;
-  pcbddc->local_mat                       = 0;
-  pcbddc->ChangeOfBasisMatrix             = 0;
-  pcbddc->usechangeofbasis                = PETSC_TRUE;
-  pcbddc->usechangeonfaces                = PETSC_FALSE;
-  pcbddc->coarse_vec                      = 0;
-  pcbddc->coarse_rhs                      = 0;
-  pcbddc->coarse_ksp                      = 0;
-  pcbddc->coarse_phi_B                    = 0;
-  pcbddc->coarse_phi_D                    = 0;
-  pcbddc->vec1_P                          = 0;
-  pcbddc->vec1_R                          = 0;
-  pcbddc->vec2_R                          = 0;
-  pcbddc->local_auxmat1                   = 0;
-  pcbddc->local_auxmat2                   = 0;
-  pcbddc->R_to_B                          = 0;
-  pcbddc->R_to_D                          = 0;
-  pcbddc->ksp_D                           = 0;
-  pcbddc->ksp_R                           = 0;
-  pcbddc->local_primal_indices            = 0;
-  pcbddc->inexact_prec_type               = PETSC_FALSE;
-  pcbddc->NeumannBoundaries               = 0;
-  pcbddc->ISForDofs                       = 0;
-  pcbddc->ISForVertices                   = 0;
-  pcbddc->n_ISForFaces                    = 0;
-  pcbddc->n_ISForEdges                    = 0;
-  pcbddc->ConstraintMatrix                = 0;
-  pcbddc->use_nnsp_true                   = PETSC_FALSE;
-  pcbddc->local_primal_sizes              = 0;
-  pcbddc->local_primal_displacements      = 0;
+  pcbddc->user_primal_vertices       = 0;
+  pcbddc->NullSpace                  = 0;
+  pcbddc->temp_solution              = 0;
+  pcbddc->original_rhs               = 0;
+  pcbddc->local_mat                  = 0;
+  pcbddc->ChangeOfBasisMatrix        = 0;
+  pcbddc->use_change_of_basis        = PETSC_TRUE;
+  pcbddc->use_change_on_faces        = PETSC_FALSE;
+  pcbddc->coarse_vec                 = 0;
+  pcbddc->coarse_rhs                 = 0;
+  pcbddc->coarse_ksp                 = 0;
+  pcbddc->coarse_phi_B               = 0;
+  pcbddc->coarse_phi_D               = 0;
+  pcbddc->vec1_P                     = 0;          
+  pcbddc->vec1_R                     = 0; 
+  pcbddc->vec2_R                     = 0; 
+  pcbddc->local_auxmat1              = 0;
+  pcbddc->local_auxmat2              = 0;
+  pcbddc->R_to_B                     = 0;
+  pcbddc->R_to_D                     = 0;
+  pcbddc->ksp_D                      = 0;
+  pcbddc->ksp_R                      = 0;
+  pcbddc->local_primal_indices       = 0;
+  pcbddc->inexact_prec_type          = PETSC_FALSE;
+  pcbddc->NeumannBoundaries          = 0;
+  pcbddc->ISForDofs                  = 0;
+  pcbddc->ConstraintMatrix           = 0;
+  pcbddc->use_nnsp_true              = PETSC_FALSE;
+  pcbddc->local_primal_sizes         = 0;
+  pcbddc->local_primal_displacements = 0;
+  pcbddc->coarse_loc_to_glob         = 0;
+  pcbddc->dbg_flag                   = PETSC_FALSE;
+  pcbddc->coarsening_ratio           = 8;
+  pcbddc->use_exact_dirichlet        = PETSC_TRUE;
+  pcbddc->current_level              = 0;
+  pcbddc->max_levels                 = 1;
   pcbddc->replicated_local_primal_indices = 0;
   pcbddc->replicated_local_primal_values  = 0;
-  pcbddc->coarse_loc_to_glob              = 0;
-  pcbddc->dbg_flag                        = PETSC_FALSE;
-  pcbddc->coarsening_ratio                = 8;
-  pcbddc->use_exact_dirichlet             = PETSC_TRUE;
-  pcbddc->current_level                   = 0;
-  pcbddc->max_levels                      = 1;
-
-  /* allocate and initialize needed graph structure */
-  ierr                      = PetscMalloc(sizeof(*mat_graph),&pcbddc->mat_graph);CHKERRQ(ierr);
-  pcbddc->mat_graph->xadj   = 0;
-  pcbddc->mat_graph->adjncy = 0;
+
+  /* create local graph structure */
+  ierr = PCBDDCGraphCreate(&pcbddc->mat_graph);CHKERRQ(ierr);  
+
+  /* scaling */
+  pcbddc->use_deluxe_scaling         = PETSC_FALSE;
+  pcbddc->work_scaling               = 0;
 
   /* function pointers */
   pc->ops->apply               = PCApply_BDDC;
   pc->ops->postsolve           = PCPostSolve_BDDC;
 
   /* composing function */
+ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",PCBDDCSetPrimalVerticesLocalIS_BDDC);CHKERRQ(ierr);
   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",PCBDDCSetCoarseningRatio_BDDC);CHKERRQ(ierr);
   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetMaxLevels_C",PCBDDCSetMaxLevels_BDDC);CHKERRQ(ierr);
   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNullSpace_C",PCBDDCSetNullSpace_BDDC);CHKERRQ(ierr);
   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseProblemType_C",PCBDDCSetCoarseProblemType_BDDC);CHKERRQ(ierr);
   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",PCBDDCSetDofsSplitting_BDDC);CHKERRQ(ierr);
   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",PCBDDCSetLocalAdjacencyGraph_BDDC);CHKERRQ(ierr);
-  ierr = PetscObjectComposeFunction((PetscObject)pc,"PCPreSolve_C",PCPreSolve_BDDC);CHKERRQ(ierr);
-  ierr = PetscObjectComposeFunction((PetscObject)pc,"PCPostSolve_C",PCPostSolve_BDDC);CHKERRQ(ierr);
   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",PCBDDCCreateFETIDPOperators_BDDC);CHKERRQ(ierr);
   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",PCBDDCMatFETIDPGetRHS_BDDC);CHKERRQ(ierr);
   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",PCBDDCMatFETIDPGetSolution_BDDC);CHKERRQ(ierr);
 /* -------------------------------------------------------------------------- */
 /* All static functions from now on                                           */
 /* -------------------------------------------------------------------------- */
-#undef __FUNCT__
-#define __FUNCT__ "PCBDDCApplyNullSpaceCorrectionPC"
-static PetscErrorCode PCBDDCApplyNullSpaceCorrectionPC(PC pc,Vec x,Vec y)
-{
-  NullSpaceCorrection_ctx *pc_ctx;
-  PetscErrorCode          ierr;
-
-  PetscFunctionBegin;
-  ierr = PCShellGetContext(pc,(void**)&pc_ctx);CHKERRQ(ierr);
-  /* E */
-  ierr = MatMultTranspose(pc_ctx->Lbasis_mat,x,pc_ctx->work_small_2);CHKERRQ(ierr);
-  ierr = MatMultAdd(pc_ctx->Kbasis_mat,pc_ctx->work_small_2,x,pc_ctx->work_full_1);CHKERRQ(ierr);
-  /* P^-1 */
-  ierr = PCApply(pc_ctx->local_pc,pc_ctx->work_full_1,pc_ctx->work_full_2);CHKERRQ(ierr);
-  /* E^T */
-  ierr = MatMultTranspose(pc_ctx->Kbasis_mat,pc_ctx->work_full_2,pc_ctx->work_small_1);CHKERRQ(ierr);
-  ierr = VecScale(pc_ctx->work_small_1,-1.0);CHKERRQ(ierr);
-  ierr = MatMultAdd(pc_ctx->Lbasis_mat,pc_ctx->work_small_1,pc_ctx->work_full_2,pc_ctx->work_full_1);CHKERRQ(ierr);
-  /* Sum contributions */
-  ierr = MatMultAdd(pc_ctx->basis_mat,pc_ctx->work_small_2,pc_ctx->work_full_1,y);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "PCBDDCDestroyNullSpaceCorrectionPC"
-static PetscErrorCode PCBDDCDestroyNullSpaceCorrectionPC(PC pc)
-{
-  NullSpaceCorrection_ctx *pc_ctx;
-  PetscErrorCode          ierr;
-
-  PetscFunctionBegin;
-  ierr = PCShellGetContext(pc,(void**)&pc_ctx);CHKERRQ(ierr);
-  ierr = VecDestroy(&pc_ctx->work_small_1);CHKERRQ(ierr);
-  ierr = VecDestroy(&pc_ctx->work_small_2);CHKERRQ(ierr);
-  ierr = VecDestroy(&pc_ctx->work_full_1);CHKERRQ(ierr);
-  ierr = VecDestroy(&pc_ctx->work_full_2);CHKERRQ(ierr);
-  ierr = MatDestroy(&pc_ctx->basis_mat);CHKERRQ(ierr);
-  ierr = MatDestroy(&pc_ctx->Lbasis_mat);CHKERRQ(ierr);
-  ierr = MatDestroy(&pc_ctx->Kbasis_mat);CHKERRQ(ierr);
-  ierr = PCDestroy(&pc_ctx->local_pc);CHKERRQ(ierr);
-  ierr = PetscFree(pc_ctx);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "PCBDDCAdaptLocalProblem"
-static PetscErrorCode PCBDDCAdaptLocalProblem(PC pc,IS local_dofs)
-{
-  extern PetscErrorCode PCBDDCApplyNullSpaceCorrectionPC(PC,Vec,Vec);
-  extern PetscErrorCode PCBDDCDestroyNullSpaceCorrectionPC(PC);
-
-  PC_BDDC                 *pcbddc = (PC_BDDC*)pc->data;
-  PC_IS                   *pcis   = (PC_IS*)pc->data;
-  Mat_IS                  * matis = (Mat_IS*)pc->pmat->data;
-  KSP                     *local_ksp;
-  PC                      newpc;
-  NullSpaceCorrection_ctx *shell_ctx;
-  Mat                     local_mat,local_pmat,small_mat,inv_small_mat;
-  MatStructure            local_mat_struct;
-  Vec                     work1,work2,work3;
-  const Vec               *nullvecs;
-  VecScatter              scatter_ctx;
-  IS                      is_aux;
-  MatFactorInfo           matinfo;
-  PetscScalar             *basis_mat,*Kbasis_mat,*array,*array_mat;
-  PetscScalar             one = 1.0,zero = 0.0, m_one = -1.0;
-  PetscInt                basis_dofs,basis_size,nnsp_size,i,k,n_I,n_R;
-  PetscBool               nnsp_has_cnst;
-  PetscErrorCode          ierr;
-
-  PetscFunctionBegin;
-  /* Infer the local solver */
-  ierr = ISGetSize(local_dofs,&basis_dofs);CHKERRQ(ierr);
-  ierr = VecGetSize(pcis->vec1_D,&n_I);CHKERRQ(ierr);
-  ierr = VecGetSize(pcbddc->vec1_R,&n_R);CHKERRQ(ierr);
-  if (basis_dofs == n_I) {
-    /* Dirichlet solver */
-    local_ksp = &pcbddc->ksp_D;
-  } else if (basis_dofs == n_R) {
-    /* Neumann solver */
-    local_ksp = &pcbddc->ksp_R;
-  } else SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in %s: unknown local IS size %d. n_I=%d, n_R=%d)\n",basis_dofs,n_I,n_R);
-  ierr = KSPGetOperators(*local_ksp,&local_mat,&local_pmat,&local_mat_struct);CHKERRQ(ierr);
-
-  /* Get null space vecs */
-  ierr = MatNullSpaceGetVecs(pcbddc->NullSpace,&nnsp_has_cnst,&nnsp_size,&nullvecs);CHKERRQ(ierr);
-  basis_size = nnsp_size;
-  if (nnsp_has_cnst) basis_size++;
-
-  /* Create shell ctx */
-  ierr = PetscMalloc(sizeof(*shell_ctx),&shell_ctx);CHKERRQ(ierr);
-
-  /* Create work vectors in shell context */
-  ierr = VecCreate(PETSC_COMM_SELF,&shell_ctx->work_small_1);CHKERRQ(ierr);
-  ierr = VecSetSizes(shell_ctx->work_small_1,basis_size,basis_size);CHKERRQ(ierr);
-  ierr = VecSetType(shell_ctx->work_small_1,VECSEQ);CHKERRQ(ierr);
-  ierr = VecDuplicate(shell_ctx->work_small_1,&shell_ctx->work_small_2);CHKERRQ(ierr);
-  ierr = VecCreate(PETSC_COMM_SELF,&shell_ctx->work_full_1);CHKERRQ(ierr);
-  ierr = VecSetSizes(shell_ctx->work_full_1,basis_dofs,basis_dofs);CHKERRQ(ierr);
-  ierr = VecSetType(shell_ctx->work_full_1,VECSEQ);CHKERRQ(ierr);
-  ierr = VecDuplicate(shell_ctx->work_full_1,&shell_ctx->work_full_2);CHKERRQ(ierr);
-
-  /* Allocate workspace */
-  ierr = MatCreateSeqDense(PETSC_COMM_SELF,basis_dofs,basis_size,NULL,&shell_ctx->basis_mat);CHKERRQ(ierr);
-  ierr = MatCreateSeqDense(PETSC_COMM_SELF,basis_dofs,basis_size,NULL,&shell_ctx->Kbasis_mat);CHKERRQ(ierr);
-  ierr = MatDenseGetArray(shell_ctx->basis_mat,&basis_mat);CHKERRQ(ierr);
-  ierr = MatDenseGetArray(shell_ctx->Kbasis_mat,&Kbasis_mat);CHKERRQ(ierr);
-
-  /* Restrict local null space on selected dofs (Dirichlet or Neumann)
-     and compute matrices N and K*N */
-  ierr = VecDuplicate(shell_ctx->work_full_1,&work1);CHKERRQ(ierr);
-  ierr = VecDuplicate(shell_ctx->work_full_1,&work2);CHKERRQ(ierr);
-  ierr = VecScatterCreate(pcis->vec1_N,local_dofs,work1,(IS)0,&scatter_ctx);CHKERRQ(ierr);
-  k    = 0;
-  for (; k<nnsp_size; k++) {
-    ierr = VecScatterBegin(matis->ctx,nullvecs[k],pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-    ierr = VecScatterEnd(matis->ctx,nullvecs[k],pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-    ierr = VecPlaceArray(work1,(const PetscScalar*)&basis_mat[k*basis_dofs]);CHKERRQ(ierr);
-    ierr = VecScatterBegin(scatter_ctx,pcis->vec1_N,work1,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-    ierr = VecScatterEnd(scatter_ctx,pcis->vec1_N,work1,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-    ierr = VecPlaceArray(work2,(const PetscScalar*)&Kbasis_mat[k*basis_dofs]);CHKERRQ(ierr);
-    ierr = MatMult(local_mat,work1,work2);CHKERRQ(ierr);
-    ierr = VecResetArray(work1);CHKERRQ(ierr);
-    ierr = VecResetArray(work2);CHKERRQ(ierr);
-  }
-  if (nnsp_has_cnst) {
-    ierr = VecPlaceArray(work1,(const PetscScalar*)&basis_mat[k*basis_dofs]);CHKERRQ(ierr);
-    ierr = VecSet(work1,one);CHKERRQ(ierr);
-    ierr = VecPlaceArray(work2,(const PetscScalar*)&Kbasis_mat[k*basis_dofs]);CHKERRQ(ierr);
-    ierr = MatMult(local_mat,work1,work2);CHKERRQ(ierr);
-    ierr = VecResetArray(work1);CHKERRQ(ierr);
-    ierr = VecResetArray(work2);CHKERRQ(ierr);
-  }
-  ierr = VecDestroy(&work1);CHKERRQ(ierr);
-  ierr = VecDestroy(&work2);CHKERRQ(ierr);
-  ierr = VecScatterDestroy(&scatter_ctx);CHKERRQ(ierr);
-  ierr = MatDenseRestoreArray(shell_ctx->basis_mat,&basis_mat);CHKERRQ(ierr);
-  ierr = MatDenseRestoreArray(shell_ctx->Kbasis_mat,&Kbasis_mat);CHKERRQ(ierr);
-
-  /* Assemble another Mat object in shell context */
-  ierr = MatTransposeMatMult(shell_ctx->basis_mat,shell_ctx->Kbasis_mat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&small_mat);CHKERRQ(ierr);
-  ierr = MatFactorInfoInitialize(&matinfo);CHKERRQ(ierr);
-  ierr = ISCreateStride(PETSC_COMM_SELF,basis_size,0,1,&is_aux);CHKERRQ(ierr);
-  ierr = MatLUFactor(small_mat,is_aux,is_aux,&matinfo);CHKERRQ(ierr);
-  ierr = ISDestroy(&is_aux);CHKERRQ(ierr);
-  ierr = PetscMalloc(basis_size*basis_size*sizeof(PetscScalar),&array_mat);CHKERRQ(ierr);
-  for (k=0; k<basis_size; k++) {
-    ierr = VecSet(shell_ctx->work_small_1,zero);CHKERRQ(ierr);
-    ierr = VecSetValue(shell_ctx->work_small_1,k,one,INSERT_VALUES);CHKERRQ(ierr);
-    ierr = VecAssemblyBegin(shell_ctx->work_small_1);CHKERRQ(ierr);
-    ierr = VecAssemblyEnd(shell_ctx->work_small_1);CHKERRQ(ierr);
-    ierr = MatSolve(small_mat,shell_ctx->work_small_1,shell_ctx->work_small_2);CHKERRQ(ierr);
-    ierr = VecGetArray(shell_ctx->work_small_2,&array);CHKERRQ(ierr);
-    for (i=0; i<basis_size; i++) array_mat[i*basis_size+k]=array[i];
-    ierr = VecRestoreArray(shell_ctx->work_small_2,&array);CHKERRQ(ierr);
-  }
-  ierr = MatCreateSeqDense(PETSC_COMM_SELF,basis_size,basis_size,array_mat,&inv_small_mat);CHKERRQ(ierr);
-  ierr = MatMatMult(shell_ctx->basis_mat,inv_small_mat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&shell_ctx->Lbasis_mat);CHKERRQ(ierr);
-  ierr = PetscFree(array_mat);CHKERRQ(ierr);
-  ierr = MatDestroy(&inv_small_mat);CHKERRQ(ierr);
-  ierr = MatDestroy(&small_mat);CHKERRQ(ierr);
-  ierr = MatScale(shell_ctx->Kbasis_mat,m_one);CHKERRQ(ierr);
-
-  /* Rebuild local PC */
-  ierr = KSPGetPC(*local_ksp,&shell_ctx->local_pc);CHKERRQ(ierr);
-  ierr = PetscObjectReference((PetscObject)shell_ctx->local_pc);CHKERRQ(ierr);
-  ierr = PCCreate(PETSC_COMM_SELF,&newpc);CHKERRQ(ierr);
-  ierr = PCSetOperators(newpc,local_mat,local_mat,SAME_PRECONDITIONER);CHKERRQ(ierr);
-  ierr = PCSetType(newpc,PCSHELL);CHKERRQ(ierr);
-  ierr = PCShellSetContext(newpc,shell_ctx);CHKERRQ(ierr);
-  ierr = PCShellSetApply(newpc,PCBDDCApplyNullSpaceCorrectionPC);CHKERRQ(ierr);
-  ierr = PCShellSetDestroy(newpc,PCBDDCDestroyNullSpaceCorrectionPC);CHKERRQ(ierr);
-  ierr = PCSetUp(newpc);CHKERRQ(ierr);
-  ierr = KSPSetPC(*local_ksp,newpc);CHKERRQ(ierr);
-  ierr = PCDestroy(&newpc);CHKERRQ(ierr);
-  ierr = KSPSetUp(*local_ksp);CHKERRQ(ierr);
-
-  /* test */
-  if (pcbddc->dbg_flag) {
-    PetscReal   test_err;
-    KSP         check_ksp;
-    PC          check_pc;
-    PetscReal   lambda_min,lambda_max;
-    Mat         test_mat;
-    PetscViewer viewer=pcbddc->dbg_viewer;
-    PetscBool   setsym,issym=PETSC_FALSE;
-
-    ierr = KSPGetPC(*local_ksp,&check_pc);CHKERRQ(ierr);
-    ierr = VecDuplicate(shell_ctx->work_full_1,&work1);CHKERRQ(ierr);
-    ierr = VecDuplicate(shell_ctx->work_full_1,&work2);CHKERRQ(ierr);
-    ierr = VecDuplicate(shell_ctx->work_full_1,&work3);CHKERRQ(ierr);
-    ierr = VecSetRandom(shell_ctx->work_small_1,NULL);CHKERRQ(ierr);
-    ierr = MatMult(shell_ctx->basis_mat,shell_ctx->work_small_1,work1);CHKERRQ(ierr);
-    ierr = VecCopy(work1,work2);CHKERRQ(ierr);
-    ierr = MatMult(local_mat,work1,work3);CHKERRQ(ierr);
-    ierr = PCApply(check_pc,work3,work1);CHKERRQ(ierr);
-    ierr = VecAXPY(work1,m_one,work2);CHKERRQ(ierr);
-    ierr = VecNorm(work1,NORM_INFINITY,&test_err);CHKERRQ(ierr);
-    ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d error for nullspace correction for ",PetscGlobalRank);
-    if (basis_dofs == n_I) {
-      ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Dirichlet ");CHKERRQ(ierr);
-    } else {
-      ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Neumann ");CHKERRQ(ierr);
-    }
-    ierr = PetscViewerASCIISynchronizedPrintf(viewer,"solver is :%1.14e\n",test_err);
-
-    ierr = MatTransposeMatMult(shell_ctx->Lbasis_mat,shell_ctx->Kbasis_mat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&test_mat);CHKERRQ(ierr);
-    ierr = MatShift(test_mat,one);CHKERRQ(ierr);
-    ierr = MatNorm(test_mat,NORM_INFINITY,&test_err);CHKERRQ(ierr);
-    ierr = MatDestroy(&test_mat);CHKERRQ(ierr);
-    ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d error for nullspace matrices is :%1.14e\n",PetscGlobalRank,test_err);
-
-    /* Create ksp object suitable for extreme eigenvalues' estimation */
-    ierr = KSPCreate(PETSC_COMM_SELF,&check_ksp);CHKERRQ(ierr);
-    ierr = KSPSetOperators(check_ksp,local_mat,local_mat,SAME_PRECONDITIONER);CHKERRQ(ierr);
-    ierr = KSPSetTolerances(check_ksp,1.e-8,1.e-8,PETSC_DEFAULT,basis_dofs);CHKERRQ(ierr);
-    ierr = KSPSetComputeSingularValues(check_ksp,PETSC_TRUE);CHKERRQ(ierr);
-    ierr = MatIsSymmetricKnown(pc->pmat,&setsym,&issym);CHKERRQ(ierr);
-    if (issym) {
-      ierr = KSPSetType(check_ksp,KSPCG);CHKERRQ(ierr);
-    }
-    ierr = KSPSetPC(check_ksp,check_pc);CHKERRQ(ierr);
-    ierr = KSPSetUp(check_ksp);CHKERRQ(ierr);
-    ierr = VecSetRandom(work1,NULL);CHKERRQ(ierr);
-    ierr = MatMult(local_mat,work1,work2);CHKERRQ(ierr);
-    ierr = KSPSolve(check_ksp,work2,work2);CHKERRQ(ierr);
-    ierr = VecAXPY(work2,m_one,work1);CHKERRQ(ierr);
-    ierr = VecNorm(work2,NORM_INFINITY,&test_err);CHKERRQ(ierr);
-    ierr = KSPComputeExtremeSingularValues(check_ksp,&lambda_max,&lambda_min);CHKERRQ(ierr);
-    ierr = KSPGetIterationNumber(check_ksp,&k);CHKERRQ(ierr);
-    ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d error for adapted KSP %1.14e (it %d, eigs %1.6e %1.6e)\n",PetscGlobalRank,test_err,k,lambda_min,lambda_max);
-    ierr = KSPDestroy(&check_ksp);CHKERRQ(ierr);
-    ierr = VecDestroy(&work1);CHKERRQ(ierr);
-    ierr = VecDestroy(&work2);CHKERRQ(ierr);
-    ierr = VecDestroy(&work3);CHKERRQ(ierr);
-  }
-  PetscFunctionReturn(0);
-}
 
-#undef __FUNCT__
+#undef __FUNCT__  
 #define __FUNCT__ "PCBDDCSetUseExactDirichlet"
 static PetscErrorCode PCBDDCSetUseExactDirichlet(PC pc,PetscBool use)
 {
-  PC_BDDC *pcbddc = (PC_BDDC*)pc->data;
+  PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;
 
   PetscFunctionBegin;
   pcbddc->use_exact_dirichlet=use;
 #define __FUNCT__ "PCBDDCSetLevel"
 static PetscErrorCode PCBDDCSetLevel(PC pc,PetscInt level)
 {
-  PC_BDDC *pcbddc = (PC_BDDC*)pc->data;
+  PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;
 
   PetscFunctionBegin;
   pcbddc->current_level=level;
   PetscFunctionReturn(0);
 }
 
-#undef __FUNCT__
-#define __FUNCT__ "PCBDDCAdaptNullSpace"
-static PetscErrorCode PCBDDCAdaptNullSpace(PC pc)
-{
-  PC_IS          *pcis   = (PC_IS*)  (pc->data);
-  PC_BDDC        *pcbddc = (PC_BDDC*)(pc->data);
-  KSP            inv_change;
-  PC             pc_change;
-  const Vec      *nsp_vecs;
-  Vec            *new_nsp_vecs;
-  PetscInt       i,nsp_size,new_nsp_size,start_new;
-  PetscBool      nsp_has_cnst;
-  MatNullSpace   new_nsp;
-  PetscErrorCode ierr;
-
-  PetscFunctionBegin;
-  ierr = MatNullSpaceGetVecs(pcbddc->NullSpace,&nsp_has_cnst,&nsp_size,&nsp_vecs);CHKERRQ(ierr);
-  ierr = KSPCreate(PETSC_COMM_SELF,&inv_change);CHKERRQ(ierr);
-  ierr = KSPSetOperators(inv_change,pcbddc->ChangeOfBasisMatrix,pcbddc->ChangeOfBasisMatrix,SAME_PRECONDITIONER);CHKERRQ(ierr);
-  ierr = KSPSetType(inv_change,KSPPREONLY);CHKERRQ(ierr);
-  ierr = KSPGetPC(inv_change,&pc_change);CHKERRQ(ierr);
-  ierr = PCSetType(pc_change,PCLU);CHKERRQ(ierr);
-  ierr = KSPSetUp(inv_change);CHKERRQ(ierr);
-
-  new_nsp_size = nsp_size;
-  if (nsp_has_cnst) new_nsp_size++;
-  ierr = PetscMalloc(new_nsp_size*sizeof(Vec),&new_nsp_vecs);CHKERRQ(ierr);
-  for (i=0;i<new_nsp_size;i++) { ierr = VecDuplicate(pcis->vec1_global,&new_nsp_vecs[i]);CHKERRQ(ierr); }
-  start_new = 0;
-  if (nsp_has_cnst) {
-    start_new = 1;
-    ierr      = VecSet(new_nsp_vecs[0],1.0);CHKERRQ(ierr);
-    ierr      = VecSet(pcis->vec1_B,1.0);CHKERRQ(ierr);
-    ierr      = KSPSolve(inv_change,pcis->vec1_B,pcis->vec1_B);
-    ierr      = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,new_nsp_vecs[0],INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
-    ierr      = VecScatterEnd  (pcis->global_to_B,pcis->vec1_B,new_nsp_vecs[0],INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
-  }
-  for (i=0; i<nsp_size; i++) {
-    ierr = VecCopy(nsp_vecs[i],new_nsp_vecs[i+start_new]);CHKERRQ(ierr);
-    ierr = VecScatterBegin(pcis->global_to_B,nsp_vecs[i],pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-    ierr = VecScatterEnd  (pcis->global_to_B,nsp_vecs[i],pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-    ierr = KSPSolve(inv_change,pcis->vec1_B,pcis->vec1_B);
-    ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,new_nsp_vecs[i+start_new],INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
-    ierr = VecScatterEnd  (pcis->global_to_B,pcis->vec1_B,new_nsp_vecs[i+start_new],INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
-  }
-  ierr = VecNormalize(new_nsp_vecs[0],NULL);CHKERRQ(ierr);
-  /* TODO : Orthonormalize vecs when new_nsp_size > 0! */
-
-  ierr = KSPDestroy(&inv_change);CHKERRQ(ierr);
-  ierr = MatNullSpaceCreate(PetscObjectComm((PetscObject)pc),PETSC_FALSE,new_nsp_size,new_nsp_vecs,&new_nsp);CHKERRQ(ierr);
-  ierr = PCBDDCSetNullSpace(pc,new_nsp);CHKERRQ(ierr);
-  ierr = MatNullSpaceDestroy(&new_nsp);CHKERRQ(ierr);
-  /*
-  ierr = MatNullSpaceTest(pcbddc->NullSpace,pc->pmat,&nsp_t);CHKERRQ(ierr);
-  printf("New Null Space, mat changed: %d\n",nsp_t);
-    temp_mat = matis->A;
-    matis->A = pcbddc->local_mat;
-    pcbddc->local_mat = temp_mat;
-  ierr = MatNullSpaceTest(pcbddc->NullSpace,pc->pmat,&nsp_t);CHKERRQ(ierr);
-  printf("New Null Space, mat original: %d\n",nsp_t);*/
-
-  for (i=0; i<new_nsp_size; i++) { ierr = VecDestroy(&new_nsp_vecs[i]);CHKERRQ(ierr); }
-  ierr = PetscFree(new_nsp_vecs);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "PCBDDCCreateFETIDPMatContext"
-static PetscErrorCode PCBDDCCreateFETIDPMatContext(PC pc, FETIDPMat_ctx **fetidpmat_ctx)
-{
-  FETIDPMat_ctx  *newctx;
-  PetscErrorCode ierr;
-
-  PetscFunctionBegin;
-  ierr = PetscMalloc(sizeof(*newctx),&newctx);CHKERRQ(ierr);
-
-  newctx->lambda_local    = 0;
-  newctx->temp_solution_B = 0;
-  newctx->temp_solution_D = 0;
-  newctx->B_delta         = 0;
-  newctx->B_Ddelta        = 0; /* theoretically belongs to the FETIDP preconditioner */
-  newctx->l2g_lambda      = 0;
-
-  /* increase the reference count for BDDC preconditioner */
-  ierr           = PetscObjectReference((PetscObject)pc);CHKERRQ(ierr);
-  newctx->pc     = pc;
-  *fetidpmat_ctx = newctx;
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "PCBDDCCreateFETIDPPCContext"
-static PetscErrorCode PCBDDCCreateFETIDPPCContext(PC pc, FETIDPPC_ctx **fetidppc_ctx)
-{
-  FETIDPPC_ctx   *newctx;
-  PetscErrorCode ierr;
-
-  PetscFunctionBegin;
-  ierr = PetscMalloc(sizeof(*newctx),&newctx);CHKERRQ(ierr);
-
-  newctx->lambda_local = 0;
-  newctx->B_Ddelta     = 0;
-  newctx->l2g_lambda   = 0;
-
-  /* increase the reference count for BDDC preconditioner */
-  ierr          = PetscObjectReference((PetscObject)pc);CHKERRQ(ierr);
-  newctx->pc    = pc;
-  *fetidppc_ctx = newctx;
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "PCBDDCDestroyFETIDPMat"
-static PetscErrorCode PCBDDCDestroyFETIDPMat(Mat A)
-{
-  FETIDPMat_ctx  *mat_ctx;
-  PetscErrorCode ierr;
-
-  PetscFunctionBegin;
-  ierr = MatShellGetContext(A,(void**)&mat_ctx);CHKERRQ(ierr);
-  ierr = VecDestroy(&mat_ctx->lambda_local);CHKERRQ(ierr);
-  ierr = VecDestroy(&mat_ctx->temp_solution_D);CHKERRQ(ierr);
-  ierr = VecDestroy(&mat_ctx->temp_solution_B);CHKERRQ(ierr);
-  ierr = MatDestroy(&mat_ctx->B_delta);CHKERRQ(ierr);
-  ierr = MatDestroy(&mat_ctx->B_Ddelta);CHKERRQ(ierr);
-  ierr = VecScatterDestroy(&mat_ctx->l2g_lambda);CHKERRQ(ierr);
-  ierr = PCDestroy(&mat_ctx->pc);CHKERRQ(ierr); /* actually it does not destroy BDDC, only decrease its reference count */
-  ierr = PetscFree(mat_ctx);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "PCBDDCDestroyFETIDPPC"
-static PetscErrorCode PCBDDCDestroyFETIDPPC(PC pc)
-{
-  FETIDPPC_ctx   *pc_ctx;
-  PetscErrorCode ierr;
+/* -------------------------------------------------------------------------- */
+#undef __FUNCT__  
+#define __FUNCT__ "PCBDDCCoarseSetUp"
+static PetscErrorCode PCBDDCCoarseSetUp(PC pc)
+{   
+  PetscErrorCode  ierr;
+
+  PC_IS*            pcis = (PC_IS*)(pc->data);
+  PC_BDDC*          pcbddc = (PC_BDDC*)pc->data;
+  Mat_IS            *matis = (Mat_IS*)pc->pmat->data;
+  Mat               change_mat_all;
+  IS                is_R_local;
+  IS                is_V_local;
+  IS                is_C_local;
+  IS                is_aux1;
+  IS                is_aux2;
+  VecType           impVecType;
+  MatType           impMatType;
+  PetscInt          n_R=0;
+  PetscInt          n_D=0;
+  PetscInt          n_B=0;
+  PetscScalar       zero=0.0;
+  PetscScalar       one=1.0;
+  PetscScalar       m_one=-1.0;
+  PetscScalar*      array;
+  PetscScalar       *coarse_submat_vals;
+  PetscInt          *idx_R_local;
+  PetscInt          *idx_V_B;
+  PetscScalar       *coarsefunctions_errors;
+  PetscScalar       *constraints_errors;
+  /* auxiliary indices */
+  PetscInt          i,j,k;
+  /* for verbose output of bddc */
+  PetscViewer       viewer=pcbddc->dbg_viewer;
+  PetscBool         dbg_flag=pcbddc->dbg_flag;
+  /* for counting coarse dofs */
+  PetscInt          n_vertices,n_constraints;
+  PetscInt          size_of_constraint;
+  PetscInt          *row_cmat_indices;
+  PetscScalar       *row_cmat_values;
+  PetscInt          *vertices,*nnz,*is_indices,*temp_indices;
+  ISLocalToGlobalMapping BtoNmap;
 
   PetscFunctionBegin;
-  ierr = PCShellGetContext(pc,(void**)&pc_ctx);CHKERRQ(ierr);
-  ierr = VecDestroy(&pc_ctx->lambda_local);CHKERRQ(ierr);
-  ierr = MatDestroy(&pc_ctx->B_Ddelta);CHKERRQ(ierr);
-  ierr = VecScatterDestroy(&pc_ctx->l2g_lambda);CHKERRQ(ierr);
-  ierr = PCDestroy(&pc_ctx->pc);CHKERRQ(ierr); /* actually it does not destroy BDDC, only decrease its reference count */
-  ierr = PetscFree(pc_ctx);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "PCBDDCSetupFETIDPMatContext"
-static PetscErrorCode PCBDDCSetupFETIDPMatContext(FETIDPMat_ctx *fetidpmat_ctx)
-{
-  PetscErrorCode ierr;
-  PC_IS          *pcis    =(PC_IS*)fetidpmat_ctx->pc->data;
-  PC_BDDC        *pcbddc  =(PC_BDDC*)fetidpmat_ctx->pc->data;
-  PCBDDCGraph    mat_graph=pcbddc->mat_graph;
-  Mat_IS         *matis   = (Mat_IS*)fetidpmat_ctx->pc->pmat->data;
-  MPI_Comm       comm     = ((PetscObject)(fetidpmat_ctx->pc))->comm;
-
-  Mat ScalingMat;
-  Vec lambda_global;
-  IS  IS_l2g_lambda;
-
-  PetscBool   skip_node,fully_redundant;
-  PetscInt    i,j,k,s,n_boundary_dofs,n_global_lambda,n_vertices,partial_sum;
-  PetscInt    n_local_lambda,n_lambda_for_dof,dual_size,n_neg_values,n_pos_values;
-  PetscMPIInt rank,nprocs;
-  PetscScalar scalar_value;
-
-  PetscInt    *vertex_indices,*temp_indices;
-  PetscInt    *dual_dofs_boundary_indices,*aux_local_numbering_1,*aux_global_numbering;
-  PetscInt    *aux_sums,*cols_B_delta,*l2g_indices;
-  PetscScalar *array,*scaling_factors,*vals_B_delta;
-  PetscInt    *aux_local_numbering_2,*dof_sizes,*dof_displs;
-  PetscInt    first_index,old_index;
-  PetscBool   first_found = PETSC_FALSE;
-
-  /* For communication of scaling factors */
-  PetscInt    *ptrs_buffer,neigh_position;
-  PetscScalar **all_factors,*send_buffer,*recv_buffer;
-  MPI_Request *send_reqs,*recv_reqs;
-
-  /* tests */
-  Vec         test_vec;
-  PetscBool   test_fetidp;
-  PetscViewer viewer;
+  /* Set Non-overlapping dimensions */
+  n_B = pcis->n_B; n_D = pcis->n - n_B;
+  /* Set types for local objects needed by BDDC precondtioner */
+  impMatType = MATSEQDENSE;
+  impVecType = VECSEQ;
+  /* get vertex indices from constraint matrix */
+  ierr = PCBDDCGetPrimalVerticesLocalIdx(pc,&n_vertices,&vertices);CHKERRQ(ierr);
+  /* Set number of constraints */
+  n_constraints = pcbddc->local_primal_size-n_vertices;
 
-  PetscFunctionBegin;
-  ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
-  ierr = MPI_Comm_size(comm,&nprocs);CHKERRQ(ierr);
-
-  /* Default type of lagrange multipliers is non-redundant */
-  fully_redundant = PETSC_FALSE;
-  ierr            = PetscOptionsGetBool(NULL,"-fetidp_fullyredundant",&fully_redundant,NULL);CHKERRQ(ierr);
-
-  /* Evaluate local and global number of lagrange multipliers */
-  ierr            = VecSet(pcis->vec1_N,0.0);CHKERRQ(ierr);
-  n_local_lambda  = 0;
-  partial_sum     = 0;
-  n_boundary_dofs = 0;
-  s               = 0;
-  n_vertices      = 0;
-  /* Get Vertices used to define the BDDC */
-  ierr = PetscMalloc(pcbddc->local_primal_size*sizeof(*vertex_indices),&vertex_indices);CHKERRQ(ierr);
-  for (i=0; i<pcbddc->local_primal_size; i++) {
-    ierr = MatGetRow(pcbddc->ConstraintMatrix,i,&j,(const PetscInt**)&temp_indices,NULL);CHKERRQ(ierr);
-    if (j == 1) {
-      vertex_indices[n_vertices]=temp_indices[0];
-      n_vertices++;
-    }
-    ierr = MatRestoreRow(pcbddc->ConstraintMatrix,i,&j,(const PetscInt**)&temp_indices,NULL);CHKERRQ(ierr);
+  /* vertices in boundary numbering */
+  ierr = PetscMalloc(n_vertices*sizeof(PetscInt),&idx_V_B);CHKERRQ(ierr);
+  ierr = ISLocalToGlobalMappingCreateIS(pcis->is_B_local,&BtoNmap);CHKERRQ(ierr);
+  ierr = ISGlobalToLocalMappingApply(BtoNmap,IS_GTOLM_DROP,n_vertices,vertices,&i,idx_V_B);CHKERRQ(ierr);
+  if (i != n_vertices) {
+    SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Error in boundary numbering for BDDC vertices! %d != %d\n",n_vertices,i);
   }
-  dual_size = pcis->n_B-n_vertices;
-
-  ierr = PetscSortInt(n_vertices,vertex_indices);CHKERRQ(ierr);
-  ierr = PetscMalloc(dual_size*sizeof(*dual_dofs_boundary_indices),&dual_dofs_boundary_indices);CHKERRQ(ierr);
-  ierr = PetscMalloc(dual_size*sizeof(*aux_local_numbering_1),&aux_local_numbering_1);CHKERRQ(ierr);
-  ierr = PetscMalloc(dual_size*sizeof(*aux_local_numbering_2),&aux_local_numbering_2);CHKERRQ(ierr);
 
-  ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr);
-  for (i=0; i<pcis->n; i++) {
-    j = mat_graph->count[i]; /* RECALL: mat_graph->count[i] does not count myself */
-    k = 0;
-    if (j > 0) k = (mat_graph->neighbours_set[i][0] == -1 ?  1 : 0);
-    j = j - k;
-    if (j > 0) n_boundary_dofs++;
-
-    skip_node = PETSC_FALSE;
-    if (s < n_vertices && vertex_indices[s]==i) { /* it works for a sorted set of vertices */
-      skip_node = PETSC_TRUE;
-      s++;
+  /* transform local matrices if needed */
+  if (pcbddc->use_change_of_basis) {
+    ierr = PetscMalloc(pcis->n*sizeof(PetscInt),&nnz);CHKERRQ(ierr);
+    ierr = ISGetIndices(pcis->is_I_local,(const PetscInt**)&is_indices);CHKERRQ(ierr);
+    for (i=0;i<n_D;i++) nnz[is_indices[i]] = 1;
+    ierr = ISRestoreIndices(pcis->is_I_local,(const PetscInt**)&is_indices);CHKERRQ(ierr);
+    ierr = ISGetIndices(pcis->is_B_local,(const PetscInt**)&is_indices);CHKERRQ(ierr);
+    k=1;
+    for (i=0;i<n_B;i++) {
+      ierr = MatGetRow(pcbddc->ChangeOfBasisMatrix,i,&j,NULL,NULL);CHKERRQ(ierr);
+      nnz[is_indices[i]]=j;
+      if (k < j) k = j;
+      ierr = MatRestoreRow(pcbddc->ChangeOfBasisMatrix,i,&j,NULL,NULL);CHKERRQ(ierr);
     }
-    if (j < 1) skip_node = PETSC_TRUE;
-    if (!skip_node) {
-      if (fully_redundant) {
-        /* fully redundant set of lagrange multipliers */
-        n_lambda_for_dof = (j*(j+1))/2;
-      } else {
-        n_lambda_for_dof = j;
-      }
-      n_local_lambda += j;
-      /* needed to evaluate global number of lagrange multipliers */
-      array[i]=(1.0*n_lambda_for_dof)/(j+1.0); /* already scaled for the next global sum */
-      /* store some data needed */
-      dual_dofs_boundary_indices[partial_sum] = n_boundary_dofs-1;
-      aux_local_numbering_1[partial_sum]      = i;
-      aux_local_numbering_2[partial_sum]      = n_lambda_for_dof;
-      partial_sum++;
+    ierr = ISRestoreIndices(pcis->is_B_local,(const PetscInt**)&is_indices);CHKERRQ(ierr);
+    /* assemble change of basis matrix on the whole set of local dofs */
+    ierr = PetscMalloc(k*sizeof(PetscInt),&temp_indices);CHKERRQ(ierr);
+    ierr = MatCreate(PETSC_COMM_SELF,&change_mat_all);CHKERRQ(ierr);
+    ierr = MatSetSizes(change_mat_all,pcis->n,pcis->n,pcis->n,pcis->n);CHKERRQ(ierr);
+    ierr = MatSetType(change_mat_all,MATSEQAIJ);CHKERRQ(ierr);
+    ierr = MatSeqAIJSetPreallocation(change_mat_all,0,nnz);CHKERRQ(ierr);
+    ierr = ISGetIndices(pcis->is_I_local,(const PetscInt**)&is_indices);CHKERRQ(ierr);
+    for (i=0;i<n_D;i++) {
+      ierr = MatSetValue(change_mat_all,is_indices[i],is_indices[i],1.0,INSERT_VALUES);CHKERRQ(ierr);
     }
-  }
-  ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr);
-
-  ierr = VecSet(pcis->vec1_global,0.0);CHKERRQ(ierr);
-  ierr = VecScatterBegin(matis->ctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
-  ierr = VecScatterEnd  (matis->ctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
-  ierr = VecSum(pcis->vec1_global,&scalar_value);CHKERRQ(ierr);
-
-  fetidpmat_ctx->n_lambda = (PetscInt) scalar_value;
-  /* printf("I found %d global multipliers (%f)\n",fetidpmat_ctx->n_lambda,scalar_value); */
-
-  /* compute global ordering of lagrange multipliers and associate l2g map */
-  ierr = VecSet(pcis->vec1_global,0.0);CHKERRQ(ierr);
-  ierr = VecSet(pcis->vec1_N,0.0);CHKERRQ(ierr);
-  ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr);
-  for (i=0;i<dual_size;i++) array[aux_local_numbering_1[i]] = aux_local_numbering_2[i];
-  ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr);
-  ierr = VecScatterBegin(matis->ctx,pcis->vec1_N,pcis->vec1_global,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
-  ierr = VecScatterEnd  (matis->ctx,pcis->vec1_N,pcis->vec1_global,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
-  ierr = VecSum(pcis->vec1_global,&scalar_value);CHKERRQ(ierr);
-  if (pcbddc->dbg_flag && (PetscInt)scalar_value != fetidpmat_ctx->n_lambda) {
-    SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Global number of multipliers mismatch! (%d!=%d)\n",(PetscInt)scalar_value,fetidpmat_ctx->n_lambda);
-  }
-
-  /* Fill pcis->vec1_global with cumulative function for global numbering */
-  ierr        = VecGetArray(pcis->vec1_global,&array);CHKERRQ(ierr);
-  ierr        = VecGetLocalSize(pcis->vec1_global,&s);CHKERRQ(ierr);
-  k           = 0;
-  first_index = -1;
-  for (i=0; i<s; i++) {
-    if (!first_found && array[i] > 0.0) {
-      first_found = PETSC_TRUE;
-      first_index = i;
-    }
-    k += (PetscInt)array[i];
-  }
-  j    = (!rank ? nprocs : 0);
-  ierr = PetscMalloc(j*sizeof(*dof_sizes),&dof_sizes);CHKERRQ(ierr);
-  ierr = PetscMalloc(j*sizeof(*dof_displs),&dof_displs);CHKERRQ(ierr);
-  ierr = MPI_Gather(&k,1,MPIU_INT,dof_sizes,1,MPIU_INT,0,comm);CHKERRQ(ierr);
-  if (!rank) {
-    dof_displs[0]=0;
-    for (i=1; i<nprocs; i++) dof_displs[i] = dof_displs[i-1]+dof_sizes[i-1];
-  }
-  ierr = MPI_Scatter(dof_displs,1,MPIU_INT,&k,1,MPIU_INT,0,comm);CHKERRQ(ierr);
-  if (first_found) {
-    array[first_index] += k;
-
-    old_index = first_index;
-    for (i=first_index+1; i<s; i++) {
-      if (array[i] > 0.0) {
-        array[i] += array[old_index];
-        old_index = i;
-      }
-    }
-  }
-  ierr = VecRestoreArray(pcis->vec1_global,&array);CHKERRQ(ierr);
-  ierr = VecSet(pcis->vec1_N,0.0);CHKERRQ(ierr);
-  ierr = VecScatterBegin(matis->ctx,pcis->vec1_global,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-  ierr = VecScatterEnd  (matis->ctx,pcis->vec1_global,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-  ierr = PetscMalloc(dual_size*sizeof(*aux_global_numbering),&aux_global_numbering);CHKERRQ(ierr);
-  ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr);
-  for (i=0; i<dual_size; i++) {
-    aux_global_numbering[i] = (PetscInt)array[aux_local_numbering_1[i]]-aux_local_numbering_2[i];
-  }
-  ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr);
-  ierr = PetscFree(aux_local_numbering_2);CHKERRQ(ierr);
-  ierr = PetscFree(dof_displs);CHKERRQ(ierr);
-  ierr = PetscFree(dof_sizes);CHKERRQ(ierr);
-
-  /* init data for scaling factors exchange */
-  partial_sum = 0;
-  j           = 0;
-
-  ierr = PetscMalloc(pcis->n_neigh*sizeof(PetscInt),&ptrs_buffer);CHKERRQ(ierr);
-  ierr = PetscMalloc((pcis->n_neigh-1)*sizeof(MPI_Request),&send_reqs);CHKERRQ(ierr);
-  ierr = PetscMalloc((pcis->n_neigh-1)*sizeof(MPI_Request),&recv_reqs);CHKERRQ(ierr);
-  ierr = PetscMalloc(pcis->n*sizeof(PetscScalar*),&all_factors);CHKERRQ(ierr);
-
-  ptrs_buffer[0] = 0;
-  for (i=1; i<pcis->n_neigh; i++) {
-    partial_sum += pcis->n_shared[i];
-    ptrs_buffer[i] = ptrs_buffer[i-1]+pcis->n_shared[i];
-  }
-  ierr = PetscMalloc(partial_sum*sizeof(PetscScalar),&send_buffer);CHKERRQ(ierr);
-  ierr = PetscMalloc(partial_sum*sizeof(PetscScalar),&recv_buffer);CHKERRQ(ierr);
-  ierr = PetscMalloc(partial_sum*sizeof(PetscScalar),&all_factors[0]);CHKERRQ(ierr);
-  for (i=0; i<pcis->n-1; i++) {
-    j = mat_graph->count[i];
-    if (j>0) {
-      k = (mat_graph->neighbours_set[i][0] == -1 ?  1 : 0);
-      j = j - k;
-    }
-    all_factors[i+1]=all_factors[i]+j;
-  }
-  /* scatter B scaling to N vec */
-  ierr = VecScatterBegin(pcis->N_to_B,pcis->D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
-  ierr = VecScatterEnd  (pcis->N_to_B,pcis->D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
-  /* communications */
-  ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr);
-  for (i=1; i<pcis->n_neigh; i++) {
-    for (j=0; j<pcis->n_shared[i]; j++)  send_buffer[ptrs_buffer[i-1]+j] = array[pcis->shared[i][j]];
-
-    j    = ptrs_buffer[i]-ptrs_buffer[i-1];
-    ierr = MPI_Isend(&send_buffer[ptrs_buffer[i-1]],j,MPIU_SCALAR,pcis->neigh[i],0,comm,&send_reqs[i-1]);CHKERRQ(ierr);
-    ierr = MPI_Irecv(&recv_buffer[ptrs_buffer[i-1]],j,MPIU_SCALAR,pcis->neigh[i],0,comm,&recv_reqs[i-1]);CHKERRQ(ierr);
-  }
-  ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr);
-  ierr = MPI_Waitall((pcis->n_neigh-1),recv_reqs,MPI_STATUSES_IGNORE);CHKERRQ(ierr);
-  /* put values in correct places */
-  for (i=1; i<pcis->n_neigh; i++) {
-    for (j=0; j<pcis->n_shared[i]; j++) {
-      k = pcis->shared[i][j];
-
-      neigh_position = 0;
-      while (mat_graph->neighbours_set[k][neigh_position] != pcis->neigh[i]) neigh_position++;
-      s = (mat_graph->neighbours_set[k][0] == -1 ? 1 : 0);
-
-      neigh_position = neigh_position - s;
-
-      all_factors[k][neigh_position]=recv_buffer[ptrs_buffer[i-1]+j];
-    }
-  }
-  ierr = MPI_Waitall((pcis->n_neigh-1),send_reqs,MPI_STATUSES_IGNORE);CHKERRQ(ierr);
-  ierr = PetscFree(send_reqs);CHKERRQ(ierr);
-  ierr = PetscFree(recv_reqs);CHKERRQ(ierr);
-  ierr = PetscFree(send_buffer);CHKERRQ(ierr);
-  ierr = PetscFree(recv_buffer);CHKERRQ(ierr);
-  ierr = PetscFree(ptrs_buffer);CHKERRQ(ierr);
-
-  /* Compute B and B_delta (local actions) */
-  ierr = PetscMalloc(pcis->n_neigh*sizeof(*aux_sums),&aux_sums);CHKERRQ(ierr);
-  ierr = PetscMalloc(n_local_lambda*sizeof(*l2g_indices),&l2g_indices);CHKERRQ(ierr);
-  ierr = PetscMalloc(n_local_lambda*sizeof(*vals_B_delta),&vals_B_delta);CHKERRQ(ierr);
-  ierr = PetscMalloc(n_local_lambda*sizeof(*cols_B_delta),&cols_B_delta);CHKERRQ(ierr);
-  ierr = PetscMalloc(n_local_lambda*sizeof(*scaling_factors),&scaling_factors);CHKERRQ(ierr);
-
-  n_global_lambda = 0;
-  partial_sum     = 0;
-
-  for (i=0;i<dual_size;i++) {
-    n_global_lambda = aux_global_numbering[i];
-    j               = mat_graph->count[aux_local_numbering_1[i]];
-    k               = (mat_graph->neighbours_set[aux_local_numbering_1[i]][0] == -1 ?  1 : 0);
-    j               = j - k;
-    aux_sums[0]     = 0;
-    for (s=1; s<j; s++) aux_sums[s]=aux_sums[s-1]+j-s+1;
-
-    array        = all_factors[aux_local_numbering_1[i]];
-    n_neg_values = 0;
-
-    while (n_neg_values < j && mat_graph->neighbours_set[aux_local_numbering_1[i]][n_neg_values+k] < rank) n_neg_values++;
-    n_pos_values = j - n_neg_values;
-
-    if (fully_redundant) {
-      for (s=0; s<n_neg_values; s++) {
-        l2g_indices    [partial_sum+s]=aux_sums[s]+n_neg_values-s-1+n_global_lambda;
-        cols_B_delta   [partial_sum+s]=dual_dofs_boundary_indices[i];
-        vals_B_delta   [partial_sum+s]=-1.0;
-        scaling_factors[partial_sum+s]=array[s];
-      }
-      for (s=0; s<n_pos_values; s++) {
-        l2g_indices    [partial_sum+s+n_neg_values]=aux_sums[n_neg_values]+s+n_global_lambda;
-        cols_B_delta   [partial_sum+s+n_neg_values]=dual_dofs_boundary_indices[i];
-        vals_B_delta   [partial_sum+s+n_neg_values]=1.0;
-        scaling_factors[partial_sum+s+n_neg_values]=array[s+n_neg_values];
-      }
-      partial_sum += j;
-    } else {
-      /* l2g_indices and default cols and vals of B_delta */
-      for (s=0; s<j; s++) {
-        l2g_indices    [partial_sum+s]=n_global_lambda+s;
-        cols_B_delta   [partial_sum+s]=dual_dofs_boundary_indices[i];
-        vals_B_delta   [partial_sum+s]=0.0;
-      }
-      /* B_delta */
-      if (n_neg_values > 0) vals_B_delta[partial_sum+n_neg_values-1] = -1.0; /* there's a rank next to me to the left */
-      if (n_neg_values < j) vals_B_delta[partial_sum+n_neg_values] = 1.0; /* there's a rank next to me to the right */
-
-      /* scaling as in Klawonn-Widlund 1999*/
-      for (s=0;s<n_neg_values;s++) {
-        scalar_value = 0.0;
-        for (k=0;k<s+1;k++) scalar_value += array[k];
-        scaling_factors[partial_sum+s] = -scalar_value;
-      }
-      for (s=0;s<n_pos_values;s++) {
-        scalar_value = 0.0;
-        for (k=s+n_neg_values;k<j;k++) scalar_value += array[k];
-        scaling_factors[partial_sum+s+n_neg_values] = scalar_value;
-      }
-      partial_sum += j;
-    }
-  }
-  ierr = PetscFree(aux_global_numbering);CHKERRQ(ierr);
-  ierr = PetscFree(aux_sums);CHKERRQ(ierr);
-  ierr = PetscFree(aux_local_numbering_1);CHKERRQ(ierr);
-  ierr = PetscFree(dual_dofs_boundary_indices);CHKERRQ(ierr);
-  ierr = PetscFree(all_factors[0]);CHKERRQ(ierr);
-  ierr = PetscFree(all_factors);CHKERRQ(ierr);
-  /* printf("I found %d local lambda dofs when numbering them (should be %d)\n",partial_sum,n_local_lambda); */
-
-  /* Local to global mapping of fetidpmat */
-  ierr = VecCreate(PETSC_COMM_SELF,&fetidpmat_ctx->lambda_local);CHKERRQ(ierr);
-  ierr = VecSetSizes(fetidpmat_ctx->lambda_local,n_local_lambda,n_local_lambda);CHKERRQ(ierr);
-  ierr = VecSetType(fetidpmat_ctx->lambda_local,VECSEQ);CHKERRQ(ierr);
-  ierr = VecCreate(comm,&lambda_global);CHKERRQ(ierr);
-  ierr = VecSetSizes(lambda_global,PETSC_DECIDE,fetidpmat_ctx->n_lambda);CHKERRQ(ierr);
-  ierr = VecSetType(lambda_global,VECMPI);CHKERRQ(ierr);
-  ierr = ISCreateGeneral(comm,n_local_lambda,l2g_indices,PETSC_OWN_POINTER,&IS_l2g_lambda);CHKERRQ(ierr);
-  ierr = VecScatterCreate(fetidpmat_ctx->lambda_local,(IS)0,lambda_global,IS_l2g_lambda,&fetidpmat_ctx->l2g_lambda);CHKERRQ(ierr);
-  ierr = ISDestroy(&IS_l2g_lambda);CHKERRQ(ierr);
-
-  /* Create local part of B_delta */
-  ierr = MatCreate(PETSC_COMM_SELF,&fetidpmat_ctx->B_delta);
-  ierr = MatSetSizes(fetidpmat_ctx->B_delta,n_local_lambda,pcis->n_B,n_local_lambda,pcis->n_B);CHKERRQ(ierr);
-  ierr = MatSetType(fetidpmat_ctx->B_delta,MATSEQAIJ);CHKERRQ(ierr);
-  ierr = MatSeqAIJSetPreallocation(fetidpmat_ctx->B_delta,1,NULL);CHKERRQ(ierr);
-  ierr = MatSetOption(fetidpmat_ctx->B_delta,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
-  for (i=0; i<n_local_lambda; i++) {
-    ierr = MatSetValue(fetidpmat_ctx->B_delta,i,cols_B_delta[i],vals_B_delta[i],INSERT_VALUES);CHKERRQ(ierr);
-  }
-  ierr = PetscFree(vals_B_delta);CHKERRQ(ierr);
-  ierr = MatAssemblyBegin(fetidpmat_ctx->B_delta,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
-  ierr = MatAssemblyEnd  (fetidpmat_ctx->B_delta,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
-
-  if (fully_redundant) {
-    ierr = MatCreate(PETSC_COMM_SELF,&ScalingMat);
-    ierr = MatSetSizes(ScalingMat,n_local_lambda,n_local_lambda,n_local_lambda,n_local_lambda);CHKERRQ(ierr);
-    ierr = MatSetType(ScalingMat,MATSEQAIJ);CHKERRQ(ierr);
-    ierr = MatSeqAIJSetPreallocation(ScalingMat,1,NULL);CHKERRQ(ierr);
-    for (i=0; i<n_local_lambda; i++) {
-      ierr = MatSetValue(ScalingMat,i,i,scaling_factors[i],INSERT_VALUES);CHKERRQ(ierr);
-    }
-    ierr = MatAssemblyBegin(ScalingMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
-    ierr = MatAssemblyEnd  (ScalingMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
-    ierr = MatMatMult(ScalingMat,fetidpmat_ctx->B_delta,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&fetidpmat_ctx->B_Ddelta);CHKERRQ(ierr);
-    ierr = MatDestroy(&ScalingMat);CHKERRQ(ierr);
-  } else {
-    ierr = MatCreate(PETSC_COMM_SELF,&fetidpmat_ctx->B_Ddelta);
-    ierr = MatSetSizes(fetidpmat_ctx->B_Ddelta,n_local_lambda,pcis->n_B,n_local_lambda,pcis->n_B);CHKERRQ(ierr);
-    ierr = MatSetType(fetidpmat_ctx->B_Ddelta,MATSEQAIJ);CHKERRQ(ierr);
-    ierr = MatSeqAIJSetPreallocation(fetidpmat_ctx->B_Ddelta,1,NULL);CHKERRQ(ierr);
-    for (i=0; i<n_local_lambda; i++) {
-      ierr = MatSetValue(fetidpmat_ctx->B_Ddelta,i,cols_B_delta[i],scaling_factors[i],INSERT_VALUES);CHKERRQ(ierr);
-    }
-    ierr = MatAssemblyBegin(fetidpmat_ctx->B_Ddelta,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
-    ierr = MatAssemblyEnd  (fetidpmat_ctx->B_Ddelta,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
-  }
-  ierr = PetscFree(scaling_factors);CHKERRQ(ierr);
-  ierr = PetscFree(cols_B_delta);CHKERRQ(ierr);
-
-  /* Create some vectors needed by fetidp */
-  ierr = VecDuplicate(pcis->vec1_B,&fetidpmat_ctx->temp_solution_B);CHKERRQ(ierr);
-  ierr = VecDuplicate(pcis->vec1_D,&fetidpmat_ctx->temp_solution_D);CHKERRQ(ierr);
-
-  test_fetidp = PETSC_FALSE;
-
-  ierr = PetscOptionsGetBool(NULL,"-fetidp_check",&test_fetidp,NULL);CHKERRQ(ierr);
-
-  if (test_fetidp) {
-
-    ierr = PetscViewerASCIIGetStdout(((PetscObject)(fetidpmat_ctx->pc))->comm,&viewer);CHKERRQ(ierr);
-    ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);CHKERRQ(ierr);
-    ierr = PetscViewerASCIIPrintf(viewer,"----------FETI_DP TESTS--------------\n");CHKERRQ(ierr);
-    ierr = PetscViewerASCIIPrintf(viewer,"All tests should return zero!\n");CHKERRQ(ierr);
-    ierr = PetscViewerASCIIPrintf(viewer,"FETIDP MAT context in the ");CHKERRQ(ierr);
-    if (fully_redundant) {
-      ierr = PetscViewerASCIIPrintf(viewer,"fully redundant case for lagrange multipliers.\n");CHKERRQ(ierr);
-    } else {
-      ierr = PetscViewerASCIIPrintf(viewer,"Non-fully redundant case for lagrange multiplier.\n");CHKERRQ(ierr);
-    }
-    ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
-
-    /* TEST A/B: Test numbering of global lambda dofs             */
-
-    ierr         = VecDuplicate(fetidpmat_ctx->lambda_local,&test_vec);CHKERRQ(ierr);
-    ierr         = VecSet(lambda_global,1.0);CHKERRQ(ierr);
-    ierr         = VecSet(test_vec,1.0);CHKERRQ(ierr);
-    ierr         = VecScatterBegin(fetidpmat_ctx->l2g_lambda,lambda_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
-    ierr         = VecScatterEnd  (fetidpmat_ctx->l2g_lambda,lambda_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
-    scalar_value = -1.0;
-    ierr         = VecAXPY(test_vec,scalar_value,fetidpmat_ctx->lambda_local);CHKERRQ(ierr);
-    ierr         = VecNorm(test_vec,NORM_INFINITY,&scalar_value);CHKERRQ(ierr);
-    ierr         = VecDestroy(&test_vec);CHKERRQ(ierr);
-    ierr         = PetscViewerASCIISynchronizedPrintf(viewer,"A[%04d]: CHECK glob to loc: % 1.14e\n",rank,scalar_value);CHKERRQ(ierr);
-    ierr         = PetscViewerFlush(viewer);CHKERRQ(ierr);
-    if (fully_redundant) {
-      ierr = VecSet(lambda_global,0.0);CHKERRQ(ierr);
-      ierr = VecSet(fetidpmat_ctx->lambda_local,0.5);CHKERRQ(ierr);
-      ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,lambda_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-      ierr = VecScatterEnd  (fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,lambda_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-      ierr = VecSum(lambda_global,&scalar_value);CHKERRQ(ierr);
-      ierr = PetscViewerASCIISynchronizedPrintf(viewer,"B[%04d]: CHECK loc to glob: % 1.14e\n",rank,scalar_value-fetidpmat_ctx->n_lambda);CHKERRQ(ierr);
-      ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
-    }
-
-    /* TEST C: It should holds B_delta*w=0, w\in\widehat{W}           */
-    /* This is the meaning of the B matrix                            */
-
-    ierr = VecSetRandom(pcis->vec1_N,NULL);CHKERRQ(ierr);
-    ierr = VecSet(pcis->vec1_global,0.0);CHKERRQ(ierr);
-    ierr = VecScatterBegin(matis->ctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
-    ierr = VecScatterEnd  (matis->ctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
-    ierr = VecScatterBegin(matis->ctx,pcis->vec1_global,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-    ierr = VecScatterEnd  (matis->ctx,pcis->vec1_global,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-    ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-    ierr = VecScatterEnd  (pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-    /* Action of B_delta */
-    ierr = MatMult(fetidpmat_ctx->B_delta,pcis->vec1_B,fetidpmat_ctx->lambda_local);CHKERRQ(ierr);
-    ierr = VecSet(lambda_global,0.0);CHKERRQ(ierr);
-    ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,lambda_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-    ierr = VecScatterEnd  (fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,lambda_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-    ierr = VecNorm(lambda_global,NORM_INFINITY,&scalar_value);CHKERRQ(ierr);
-    ierr = PetscViewerASCIIPrintf(viewer,"C[coll]: CHECK infty norm of B_delta*w (w continuous): % 1.14e\n",scalar_value);CHKERRQ(ierr);
-    ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
-
-    /* TEST D: It should holds E_Dw = w - P_Dw w\in\widetilde{W}     */
-    /* E_D = R_D^TR                                                   */
-    /* P_D = B_{D,delta}^T B_{delta}                                  */
-    /* eq.44 Mandel Tezaur and Dohrmann 2005                          */
-
-    /* compute a random vector in \widetilde{W} */
-    ierr = VecSetRandom(pcis->vec1_N,NULL);CHKERRQ(ierr);
-
-    scalar_value = 0.0; /* set zero at vertices */
-    ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr);
-    for (i=0;i<n_vertices;i++) array[vertex_indices[i]] = scalar_value;
-    ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr);
-
-    /* store w for final comparison */
-    ierr = VecDuplicate(pcis->vec1_B,&test_vec);CHKERRQ(ierr);
-    ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_N,test_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-    ierr = VecScatterEnd  (pcis->N_to_B,pcis->vec1_N,test_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-
-    /* Jump operator P_D : results stored in pcis->vec1_B */
-
-    ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-    ierr = VecScatterEnd  (pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-    /* Action of B_delta */
-    ierr = MatMult(fetidpmat_ctx->B_delta,pcis->vec1_B,fetidpmat_ctx->lambda_local);CHKERRQ(ierr);
-    ierr = VecSet(lambda_global,0.0);CHKERRQ(ierr);