xref: /petsc/src/ksp/pc/impls/bddc/bddcschurs.c (revision 08122e4325a89b7e1c4fe454da4f11b2780b06c8) 
1 #include <../src/ksp/pc/impls/bddc/bddc.h>
2 #include <../src/ksp/pc/impls/bddc/bddcprivate.h>
3 #include <petscblaslapack.h>
4 
5 PETSC_STATIC_INLINE PetscErrorCode PCBDDCAdjGetNextLayer_Private(PetscInt*,PetscInt,PetscBT,PetscInt*,PetscInt*,PetscInt*);
6 static PetscErrorCode PCBDDCComputeExplicitSchur(Mat M, Mat *S);
7 
8 #undef __FUNCT__
9 #define __FUNCT__ "PCBDDCComputeExplicitSchur"
10 static PetscErrorCode PCBDDCComputeExplicitSchur(Mat M, Mat *S)
11 {
12   Mat            B, C, D, Bd, Cd, AinvBd;
13   KSP            ksp;
14   PC             pc;
15   PetscBool      isLU, isILU, isCHOL, Bdense, Cdense;
16   PetscReal      fill = 2.0;
17   PetscMPIInt    size;
18   PetscErrorCode ierr;
19 
20   PetscFunctionBegin;
21   ierr = MPI_Comm_size(PetscObjectComm((PetscObject)M),&size);CHKERRQ(ierr);
22   if (size != 1) {
23     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for parallel matrices");
24   }
25   ierr = MatSchurComplementGetSubMatrices(M, NULL, NULL, &B, &C, &D);CHKERRQ(ierr);
26   ierr = MatSchurComplementGetKSP(M, &ksp);CHKERRQ(ierr);
27   ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr);
28   ierr = PetscObjectTypeCompare((PetscObject) pc, PCLU, &isLU);CHKERRQ(ierr);
29   ierr = PetscObjectTypeCompare((PetscObject) pc, PCILU, &isILU);CHKERRQ(ierr);
30   ierr = PetscObjectTypeCompare((PetscObject) pc, PCCHOLESKY, &isCHOL);CHKERRQ(ierr);
31   ierr = PetscObjectTypeCompare((PetscObject) B, MATSEQDENSE, &Bdense);CHKERRQ(ierr);
32   ierr = PetscObjectTypeCompare((PetscObject) C, MATSEQDENSE, &Cdense);CHKERRQ(ierr);
33   if (!Bdense) {
34     ierr = MatConvert(B, MATSEQDENSE, MAT_INITIAL_MATRIX, &Bd);CHKERRQ(ierr);
35   } else {
36     Bd = B;
37   }
38 
39   if (isLU || isILU || isCHOL) {
40     Mat fact;
41 
42     ierr = KSPSetUp(ksp);CHKERRQ(ierr);
43     ierr = PCFactorGetMatrix(pc, &fact);CHKERRQ(ierr);
44     ierr = MatDuplicate(Bd, MAT_DO_NOT_COPY_VALUES, &AinvBd);CHKERRQ(ierr);
45     ierr = MatMatSolve(fact, Bd, AinvBd);CHKERRQ(ierr);
46   } else {
47     Mat Ainvd;
48 
49     ierr = PCComputeExplicitOperator(pc, &Ainvd);CHKERRQ(ierr);
50     ierr = MatMatMult(Ainvd, Bd, MAT_INITIAL_MATRIX, fill, &AinvBd);CHKERRQ(ierr);
51     ierr = MatDestroy(&Ainvd);CHKERRQ(ierr);
52   }
53   if (!Bdense) {
54     ierr = MatDestroy(&Bd);CHKERRQ(ierr);
55   }
56   if (!Cdense) {
57     ierr = MatConvert(C, MATSEQDENSE, MAT_INITIAL_MATRIX, &Cd);CHKERRQ(ierr);
58   } else {
59     Cd = C;
60   }
61 
62   ierr = MatMatMult(Cd, AinvBd, MAT_INITIAL_MATRIX, fill, S);CHKERRQ(ierr);
63   ierr = MatDestroy(&AinvBd);CHKERRQ(ierr);
64   if (!Cdense) {
65     ierr = MatDestroy(&Cd);CHKERRQ(ierr);
66   }
67 
68   if (D) {
69     Mat       Dd;
70     PetscBool Ddense;
71 
72     ierr = PetscObjectTypeCompare((PetscObject)D,MATSEQDENSE,&Ddense);CHKERRQ(ierr);
73     if (!Ddense) {
74       ierr = MatConvert(D, MATSEQDENSE, MAT_INITIAL_MATRIX, &Dd);CHKERRQ(ierr);
75     } else {
76       Dd = D;
77     }
78     ierr = MatAYPX(*S,-1.0,Dd,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
79     if (!Ddense) {
80       ierr = MatDestroy(&Dd);CHKERRQ(ierr);
81     }
82   } else {
83     ierr = MatScale(*S,-1.0);CHKERRQ(ierr);
84   }
85   PetscFunctionReturn(0);
86 }
87 
88 #undef __FUNCT__
89 #define __FUNCT__ "PCBDDCSubSchursSetUp"
90 PetscErrorCode PCBDDCSubSchursSetUp(PCBDDCSubSchurs sub_schurs, PetscInt xadj[], PetscInt adjncy[], PetscInt nlayers, PetscBool compute_Stilda, PetscBool deluxe, PetscBool use_edges, PetscBool use_faces)
91 {
92   Mat                    A_II,A_IB,A_BI,A_BB;
93   Mat                    AE_II,*AE_IE,*AE_EI,*AE_EE;
94   Mat                    S_all,global_schur_subsets,work_mat;
95   Mat                    S_Ej_tilda_all,S_Ej_inv_all;
96   ISLocalToGlobalMapping l2gmap_subsets;
97   IS                     is_I,*is_subset_B,temp_is;
98   PetscInt               *nnz,*all_local_idx_G,*all_local_idx_B,*all_local_idx_N;
99   PetscInt               i,subset_size,max_subset_size;
100   PetscInt               extra,local_size,global_size;
101   PetscBool              is_symmetric,Stilda_computed;
102   PetscBLASInt           B_N,B_ierr,B_lwork,*pivots;
103   PetscScalar            *work;
104   PetscErrorCode         ierr;
105 
106   PetscFunctionBegin;
107   /* get Schur complement matrices */
108   if (!sub_schurs->use_mumps) {
109     if (compute_Stilda) {
110       SETERRQ(PetscObjectComm((PetscObject)sub_schurs->l2gmap),PETSC_ERR_SUP,"Computation of Stildas requires MUMPS");
111     }
112     ierr = MatSchurComplementGetSubMatrices(sub_schurs->S,&A_II,NULL,&A_IB,&A_BI,&A_BB);CHKERRQ(ierr);
113     ierr = PetscMalloc4(sub_schurs->n_subs,&is_subset_B,
114                         sub_schurs->n_subs,&AE_IE,
115                         sub_schurs->n_subs,&AE_EI,
116                         sub_schurs->n_subs,&AE_EE);CHKERRQ(ierr);
117   }
118 
119   /* determine interior problems */
120   ierr = ISDestroy(&sub_schurs->is_I_layer);CHKERRQ(ierr);
121   if (nlayers >= 0 && xadj != NULL && adjncy != NULL) { /* Interior problems can be different from the original one */
122     PetscBT                touched;
123     const PetscInt*        idx_B;
124     PetscInt               n_I,n_B,n_local_dofs,n_prev_added,j,layer,*local_numbering;
125 
126     /* get sizes */
127     ierr = ISGetLocalSize(sub_schurs->is_I,&n_I);CHKERRQ(ierr);
128     ierr = ISGetLocalSize(sub_schurs->is_B,&n_B);CHKERRQ(ierr);
129 
130     ierr = PetscMalloc1(n_I+n_B,&local_numbering);CHKERRQ(ierr);
131     ierr = PetscBTCreate(n_I+n_B,&touched);CHKERRQ(ierr);
132     ierr = PetscBTMemzero(n_I+n_B,touched);CHKERRQ(ierr);
133 
134     /* all boundary dofs must be skipped when adding layers */
135     ierr = ISGetIndices(sub_schurs->is_B,&idx_B);CHKERRQ(ierr);
136     for (j=0;j<n_B;j++) {
137       ierr = PetscBTSet(touched,idx_B[j]);CHKERRQ(ierr);
138     }
139     ierr = PetscMemcpy(local_numbering,idx_B,n_B*sizeof(PetscInt));CHKERRQ(ierr);
140     ierr = ISRestoreIndices(sub_schurs->is_B,&idx_B);CHKERRQ(ierr);
141 
142     /* add prescribed number of layers of dofs */
143     n_local_dofs = n_B;
144     n_prev_added = n_B;
145     for (layer=0;layer<nlayers;layer++) {
146       PetscInt n_added;
147       if (n_local_dofs == n_I+n_B) break;
148       if (n_local_dofs > n_I+n_B) {
149         SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error querying layer %d. Out of bound access (%d > %d)",layer,n_local_dofs,n_I+n_B);
150       }
151       ierr = PCBDDCAdjGetNextLayer_Private(local_numbering+n_local_dofs,n_prev_added,touched,xadj,adjncy,&n_added);CHKERRQ(ierr);
152       n_prev_added = n_added;
153       n_local_dofs += n_added;
154       if (!n_added) break;
155     }
156     ierr = PetscBTDestroy(&touched);CHKERRQ(ierr);
157 
158     /* IS for I layer dofs in original numbering */
159     ierr = ISCreateGeneral(PetscObjectComm((PetscObject)sub_schurs->is_I),n_local_dofs-n_B,local_numbering+n_B,PETSC_COPY_VALUES,&sub_schurs->is_I_layer);CHKERRQ(ierr);
160     ierr = PetscFree(local_numbering);CHKERRQ(ierr);
161     ierr = ISSort(sub_schurs->is_I_layer);CHKERRQ(ierr);
162     /* IS for I layer dofs in I numbering */
163     if (!sub_schurs->use_mumps) {
164       ISLocalToGlobalMapping ItoNmap;
165       ierr = ISLocalToGlobalMappingCreateIS(sub_schurs->is_I,&ItoNmap);CHKERRQ(ierr);
166       ierr = ISGlobalToLocalMappingApplyIS(ItoNmap,IS_GTOLM_DROP,sub_schurs->is_I_layer,&is_I);CHKERRQ(ierr);
167       ierr = ISLocalToGlobalMappingDestroy(&ItoNmap);CHKERRQ(ierr);
168 
169       /* II block */
170       ierr = MatGetSubMatrix(A_II,is_I,is_I,MAT_INITIAL_MATRIX,&AE_II);CHKERRQ(ierr);
171     }
172   } else {
173     PetscInt n_I;
174 
175     /* IS for I dofs in original numbering */
176     ierr = PetscObjectReference((PetscObject)sub_schurs->is_I);CHKERRQ(ierr);
177     sub_schurs->is_I_layer = sub_schurs->is_I;
178 
179     /* IS for I dofs in I numbering (strided 1) */
180     if (!sub_schurs->use_mumps) {
181       ierr = ISGetSize(sub_schurs->is_I,&n_I);CHKERRQ(ierr);
182       ierr = ISCreateStride(PetscObjectComm((PetscObject)sub_schurs->is_I),n_I,0,1,&is_I);CHKERRQ(ierr);
183 
184       /* II block is the same */
185       ierr = PetscObjectReference((PetscObject)A_II);CHKERRQ(ierr);
186       AE_II = A_II;
187     }
188   }
189 
190   /* Get info on subset sizes and sum of all subsets sizes */
191   max_subset_size = 0;
192   local_size = 0;
193   for (i=0;i<sub_schurs->n_subs_seq;i++) {
194     PetscInt j = sub_schurs->index_sequential[i];
195     ierr = ISGetLocalSize(sub_schurs->is_subs[j],&subset_size);CHKERRQ(ierr);
196     max_subset_size = PetscMax(subset_size,max_subset_size);
197     local_size += subset_size;
198   }
199 
200   /* Work arrays for local indices */
201   ierr = PetscMalloc1(local_size,&all_local_idx_B);CHKERRQ(ierr);
202   extra = 0;
203   if (sub_schurs->use_mumps) {
204     ierr = ISGetLocalSize(sub_schurs->is_I_layer,&extra);CHKERRQ(ierr);
205   }
206   ierr = PetscMalloc1(local_size+extra,&all_local_idx_N);CHKERRQ(ierr);
207   if (extra) {
208     const PetscInt *idxs;
209     ierr = ISGetIndices(sub_schurs->is_I_layer,&idxs);CHKERRQ(ierr);
210     ierr = PetscMemcpy(all_local_idx_N,idxs,extra*sizeof(PetscInt));CHKERRQ(ierr);
211     ierr = ISRestoreIndices(sub_schurs->is_I_layer,&idxs);CHKERRQ(ierr);
212   }
213   ierr = PetscMalloc1(local_size,&nnz);CHKERRQ(ierr);
214 
215   /* Get local indices in local numbering */
216   local_size = 0;
217   for (i=0;i<sub_schurs->n_subs_seq;i++) {
218     PetscInt j;
219     const    PetscInt *idxs;
220 
221     PetscInt local_problem_index = sub_schurs->index_sequential[i];
222     ierr = ISGetLocalSize(sub_schurs->is_subs[local_problem_index],&subset_size);CHKERRQ(ierr);
223     ierr = ISGetIndices(sub_schurs->is_subs[local_problem_index],&idxs);CHKERRQ(ierr);
224     /* subset indices in local numbering */
225     ierr = PetscMemcpy(all_local_idx_N+local_size+extra,idxs,subset_size*sizeof(PetscInt));CHKERRQ(ierr);
226     ierr = ISRestoreIndices(sub_schurs->is_subs[local_problem_index],&idxs);CHKERRQ(ierr);
227     for (j=0;j<subset_size;j++) nnz[local_size+j] = subset_size;
228     local_size += subset_size;
229   }
230 
231   S_all = NULL;
232   ierr = MatIsSymmetric(sub_schurs->A,0.0,&is_symmetric);CHKERRQ(ierr);
233   if (!sub_schurs->use_mumps) {
234     /* subsets in original and boundary numbering */
235     for (i=0;i<sub_schurs->n_subs;i++) {
236       ierr = ISGlobalToLocalMappingApplyIS(sub_schurs->BtoNmap,IS_GTOLM_DROP,sub_schurs->is_subs[i],&is_subset_B[i]);CHKERRQ(ierr);
237     }
238 
239     /* EE block */
240     for (i=0;i<sub_schurs->n_subs;i++) {
241       ierr = MatGetSubMatrix(A_BB,is_subset_B[i],is_subset_B[i],MAT_INITIAL_MATRIX,&AE_EE[i]);CHKERRQ(ierr);
242     }
243     /* IE block */
244     for (i=0;i<sub_schurs->n_subs;i++) {
245       ierr = MatGetSubMatrix(A_IB,is_I,is_subset_B[i],MAT_INITIAL_MATRIX,&AE_IE[i]);CHKERRQ(ierr);
246     }
247     /* EI block */
248     for (i=0;i<sub_schurs->n_subs;i++) {
249       ierr = MatGetSubMatrix(A_BI,is_subset_B[i],is_I,MAT_INITIAL_MATRIX,&AE_EI[i]);CHKERRQ(ierr);
250     }
251 
252     /* setup Schur complements on subset */
253     for (i=0;i<sub_schurs->n_subs;i++) {
254       ierr = MatDestroy(&sub_schurs->S_Ej[i]);CHKERRQ(ierr);
255       ierr = MatCreateSchurComplement(AE_II,AE_II,AE_IE[i],AE_EI[i],AE_EE[i],&sub_schurs->S_Ej[i]);CHKERRQ(ierr);
256       if (AE_II == A_II) { /* we can reuse the same ksp */
257         KSP ksp;
258         ierr = MatSchurComplementGetKSP(sub_schurs->S,&ksp);CHKERRQ(ierr);
259         ierr = MatSchurComplementSetKSP(sub_schurs->S_Ej[i],ksp);CHKERRQ(ierr);
260       } else { /* build new ksp object which inherits ksp and pc types from the original one */
261         KSP      origksp,schurksp;
262         PC       origpc,schurpc;
263         KSPType  ksp_type;
264         PCType   pc_type;
265         PetscInt n_internal;
266 
267         ierr = MatSchurComplementGetKSP(sub_schurs->S,&origksp);CHKERRQ(ierr);
268         ierr = MatSchurComplementGetKSP(sub_schurs->S_Ej[i],&schurksp);CHKERRQ(ierr);
269         ierr = KSPGetType(origksp,&ksp_type);CHKERRQ(ierr);
270         ierr = KSPSetType(schurksp,ksp_type);CHKERRQ(ierr);
271         ierr = KSPGetPC(schurksp,&schurpc);CHKERRQ(ierr);
272         ierr = KSPGetPC(origksp,&origpc);CHKERRQ(ierr);
273         ierr = PCGetType(origpc,&pc_type);CHKERRQ(ierr);
274         ierr = PCSetType(schurpc,pc_type);CHKERRQ(ierr);
275         ierr = ISGetSize(is_I,&n_internal);CHKERRQ(ierr);
276         if (n_internal) { /* UMFPACK gives error with 0 sized problems */
277           MatSolverPackage solver=NULL;
278           ierr = PCFactorGetMatSolverPackage(origpc,(const MatSolverPackage*)&solver);CHKERRQ(ierr);
279           if (solver) {
280             ierr = PCFactorSetMatSolverPackage(schurpc,solver);CHKERRQ(ierr);
281           }
282         }
283         ierr = KSPSetUp(schurksp);CHKERRQ(ierr);
284       }
285     }
286     /* free */
287     ierr = ISDestroy(&is_I);CHKERRQ(ierr);
288     ierr = MatDestroy(&AE_II);CHKERRQ(ierr);
289     for (i=0;i<sub_schurs->n_subs;i++) {
290       ierr = MatDestroy(&AE_EE[i]);CHKERRQ(ierr);
291       ierr = MatDestroy(&AE_IE[i]);CHKERRQ(ierr);
292       ierr = MatDestroy(&AE_EI[i]);CHKERRQ(ierr);
293       ierr = ISDestroy(&is_subset_B[i]);CHKERRQ(ierr);
294     }
295     ierr = PetscFree4(is_subset_B,AE_IE,AE_EI,AE_EE);CHKERRQ(ierr);
296 #if defined(PETSC_HAVE_MUMPS)
297   } else {
298     Mat           A,F;
299     IS            is_A_all;
300     PetscInt      *idxs_schur,n_I;
301 
302     /* get working mat */
303     ierr = ISGetLocalSize(sub_schurs->is_I_layer,&n_I);CHKERRQ(ierr);
304     ierr = ISCreateGeneral(PETSC_COMM_SELF,local_size+n_I,all_local_idx_N,PETSC_COPY_VALUES,&is_A_all);CHKERRQ(ierr);
305     ierr = MatGetSubMatrixUnsorted(sub_schurs->A,is_A_all,is_A_all,MAT_INITIAL_MATRIX,&A);CHKERRQ(ierr);
306     ierr = ISDestroy(&is_A_all);CHKERRQ(ierr);
307 
308     if (n_I) {
309       if (is_symmetric) {
310         ierr = MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_CHOLESKY,&F);CHKERRQ(ierr);
311       } else {
312         ierr = MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_LU,&F);CHKERRQ(ierr);
313       }
314 
315       /* subsets ordered last */
316       ierr = PetscMalloc1(local_size,&idxs_schur);CHKERRQ(ierr);
317       for (i=0;i<local_size;i++) {
318         idxs_schur[i] = n_I+i+1;
319       }
320       ierr = MatMumpsSetSchurIndices(F,local_size,idxs_schur);CHKERRQ(ierr);
321       ierr = PetscFree(idxs_schur);CHKERRQ(ierr);
322 
323       /* factorization step */
324       if (is_symmetric) {
325         ierr = MatCholeskyFactorSymbolic(F,A,NULL,NULL);CHKERRQ(ierr);
326         ierr = MatCholeskyFactorNumeric(F,A,NULL);CHKERRQ(ierr);
327       } else {
328         ierr = MatLUFactorSymbolic(F,A,NULL,NULL,NULL);CHKERRQ(ierr);
329         ierr = MatLUFactorNumeric(F,A,NULL);CHKERRQ(ierr);
330       }
331 
332       /* get explicit Schur Complement computed during numeric factorization */
333       ierr = MatMumpsGetSchurComplement(F,&S_all);CHKERRQ(ierr);
334       ierr = MatDestroy(&F);CHKERRQ(ierr);
335     } else {
336       ierr = MatConvert(A,MATSEQDENSE,MAT_INITIAL_MATRIX,&S_all);CHKERRQ(ierr);
337     }
338     ierr = MatDestroy(&A);CHKERRQ(ierr);
339 #endif
340   }
341 
342   if (!sub_schurs->n_subs_seq_g) {
343     PetscFunctionReturn(0);
344   }
345 
346   /* subset indices in local boundary numbering */
347   if (!sub_schurs->is_Ej_all) {
348     ierr = ISGlobalToLocalMappingApply(sub_schurs->BtoNmap,IS_GTOLM_DROP,local_size,all_local_idx_N+extra,&subset_size,all_local_idx_B);CHKERRQ(ierr);
349     if (subset_size != local_size) {
350       SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in sub_schurs serial (BtoNmap)! %d != %d\n",subset_size,local_size);
351     }
352     ierr = ISCreateGeneral(PETSC_COMM_SELF,local_size,all_local_idx_B,PETSC_OWN_POINTER,&sub_schurs->is_Ej_all);CHKERRQ(ierr);
353   }
354 
355   /* Local matrix of all local Schur on subsets transposed */
356   if (!sub_schurs->S_Ej_all) {
357     ierr = MatCreate(PETSC_COMM_SELF,&sub_schurs->S_Ej_all);CHKERRQ(ierr);
358     ierr = MatSetSizes(sub_schurs->S_Ej_all,PETSC_DECIDE,PETSC_DECIDE,local_size,local_size);CHKERRQ(ierr);
359     ierr = MatSetType(sub_schurs->S_Ej_all,MATAIJ);CHKERRQ(ierr);
360     ierr = MatSeqAIJSetPreallocation(sub_schurs->S_Ej_all,0,nnz);CHKERRQ(ierr);
361   } else {
362     ierr = MatZeroEntries(sub_schurs->S_Ej_all);CHKERRQ(ierr);
363   }
364 
365   S_Ej_tilda_all = 0;
366   S_Ej_inv_all = 0;
367   if (sub_schurs->n_subs && deluxe) { /* workspace needed only for GETRI */
368     PetscScalar lwork;
369 
370     B_lwork = -1;
371     ierr = PetscBLASIntCast(max_subset_size,&B_N);CHKERRQ(ierr);
372     ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
373     PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&B_N,work,&B_N,pivots,&lwork,&B_lwork,&B_ierr));
374     ierr = PetscFPTrapPop();CHKERRQ(ierr);
375     if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in query to GETRI Lapack routine %d",(int)B_ierr);
376     ierr = PetscBLASIntCast((PetscInt)lwork,&B_lwork);CHKERRQ(ierr);
377     ierr = PetscMalloc2(B_lwork,&work,max_subset_size,&pivots);CHKERRQ(ierr);
378   } else {
379     work = NULL;
380     pivots = NULL;
381   }
382 
383   ierr = PetscBTMemzero(sub_schurs->n_subs,sub_schurs->computed_Stilda_subs);CHKERRQ(ierr);
384   if (!sub_schurs->use_mumps) {
385     PetscScalar *fill_vals;
386     PetscInt    *dummy_idx;
387 
388     /* Work arrays */
389     ierr = PetscMalloc1(max_subset_size,&dummy_idx);CHKERRQ(ierr);
390     /* Loop on local problems end compute Schur complements explicitly */
391     local_size = 0;
392     for (i=0;i<sub_schurs->n_subs_seq;i++) {
393       Mat       S_Ej_expl;
394       PetscInt  j,lpi = sub_schurs->index_sequential[i];
395       PetscBool Sdense;
396 
397       ierr = PCBDDCComputeExplicitSchur(sub_schurs->S_Ej[lpi],&S_Ej_expl);CHKERRQ(ierr);
398       ierr = ISGetLocalSize(sub_schurs->is_subs[lpi],&subset_size);CHKERRQ(ierr);
399       ierr = PetscObjectTypeCompare((PetscObject)S_Ej_expl,MATSEQDENSE,&Sdense);CHKERRQ(ierr);
400       if (Sdense) {
401         for (j=0;j<subset_size;j++) {
402           dummy_idx[j]=local_size+j;
403         }
404         ierr = MatDenseGetArray(S_Ej_expl,&fill_vals);CHKERRQ(ierr);
405         ierr = MatSetValues(sub_schurs->S_Ej_all,subset_size,dummy_idx,subset_size,dummy_idx,fill_vals,INSERT_VALUES);CHKERRQ(ierr);
406         ierr = MatDenseRestoreArray(S_Ej_expl,&fill_vals);CHKERRQ(ierr);
407       } else {
408         SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not yet implemented for sparse matrices");
409       }
410       local_size += subset_size;
411       ierr = MatDestroy(&sub_schurs->S_Ej[lpi]);CHKERRQ(ierr);
412       sub_schurs->S_Ej[lpi] = S_Ej_expl;
413     }
414     /* Stildas are not computed without mumps */
415     for (i=0;i<sub_schurs->n_subs;i++) {
416       ierr = MatDestroy(&sub_schurs->S_Ej_tilda[i]);CHKERRQ(ierr);
417     }
418     ierr = PetscFree(dummy_idx);CHKERRQ(ierr);
419   } else {
420     PetscInt  *dummy_idx,n_all;
421 
422     if (compute_Stilda) {
423       ierr = MatCreate(PETSC_COMM_SELF,&S_Ej_tilda_all);CHKERRQ(ierr);
424       ierr = MatSetSizes(S_Ej_tilda_all,PETSC_DECIDE,PETSC_DECIDE,local_size,local_size);CHKERRQ(ierr);
425       ierr = MatSetType(S_Ej_tilda_all,MATAIJ);CHKERRQ(ierr);
426       ierr = MatSeqAIJSetPreallocation(S_Ej_tilda_all,0,nnz);CHKERRQ(ierr);
427       if (deluxe) {
428         ierr = MatCreate(PETSC_COMM_SELF,&S_Ej_inv_all);CHKERRQ(ierr);
429         ierr = MatSetSizes(S_Ej_inv_all,PETSC_DECIDE,PETSC_DECIDE,local_size,local_size);CHKERRQ(ierr);
430         ierr = MatSetType(S_Ej_inv_all,MATAIJ);CHKERRQ(ierr);
431         ierr = MatSeqAIJSetPreallocation(S_Ej_inv_all,0,nnz);CHKERRQ(ierr);
432       }
433     }
434 
435     ierr = MatGetSize(S_all,&n_all,NULL);CHKERRQ(ierr);
436     local_size = 0;
437     /* Work arrays */
438     ierr = PetscMalloc1(max_subset_size,&dummy_idx);CHKERRQ(ierr);
439 
440     Stilda_computed = PETSC_FALSE;
441     for (i=0;i<sub_schurs->n_subs;i++) {
442       IS  is_E;
443       PetscScalar *vals;
444       PetscInt j;
445 
446       ierr = ISGetLocalSize(sub_schurs->is_subs[i],&subset_size);CHKERRQ(ierr);
447       for (j=0;j<subset_size;j++) {
448         dummy_idx[j]=local_size+j;
449       }
450       ierr = ISCreateStride(PetscObjectComm((PetscObject)sub_schurs->is_I),subset_size,local_size,1,&is_E);CHKERRQ(ierr);
451       if (sub_schurs->S_Ej[i]) {
452         ierr = MatGetSubMatrix(S_all,is_E,is_E,MAT_REUSE_MATRIX,&sub_schurs->S_Ej[i]);CHKERRQ(ierr);
453       } else {
454         ierr = MatGetSubMatrix(S_all,is_E,is_E,MAT_INITIAL_MATRIX,&sub_schurs->S_Ej[i]);CHKERRQ(ierr);
455       }
456       ierr = MatDenseGetArray(sub_schurs->S_Ej[i],&vals);CHKERRQ(ierr);
457       ierr = MatSetValues(sub_schurs->S_Ej_all,subset_size,dummy_idx,subset_size,dummy_idx,vals,INSERT_VALUES);CHKERRQ(ierr);
458       ierr = MatDenseRestoreArray(sub_schurs->S_Ej[i],&vals);CHKERRQ(ierr);
459 
460       if ((PetscBTLookup(sub_schurs->is_edge,i) && use_edges) || use_faces) {
461         ierr = PetscBTSet(sub_schurs->computed_Stilda_subs,i);CHKERRQ(ierr);
462         Stilda_computed = PETSC_TRUE;
463       }
464 
465       if (PetscBTLookup(sub_schurs->computed_Stilda_subs,i)) {
466         Mat Stilda;
467         if (sub_schurs->n_subs == 1) {
468           ierr = PetscObjectReference((PetscObject)sub_schurs->S_Ej[i]);CHKERRQ(ierr);
469           Stilda = sub_schurs->S_Ej[i];
470         } else {
471           KSP ksp;
472           PC  pc;
473           Mat S_EF,S_FE,S_FF,Stilda_impl;
474           IS  is_F;
475           PetscScalar eps=1.e-8;
476           PetscBool chop=PETSC_FALSE;
477 
478           ierr = ISComplement(is_E,0,n_all,&is_F);CHKERRQ(ierr);
479           ierr = MatGetSubMatrix(S_all,is_E,is_F,MAT_INITIAL_MATRIX,&S_EF);CHKERRQ(ierr);
480           ierr = MatGetSubMatrix(S_all,is_F,is_F,MAT_INITIAL_MATRIX,&S_FF);CHKERRQ(ierr);
481           ierr = MatGetSubMatrix(S_all,is_F,is_E,MAT_INITIAL_MATRIX,&S_FE);CHKERRQ(ierr);
482           ierr = ISDestroy(&is_F);CHKERRQ(ierr);
483           if (chop) {
484             ierr = MatChop(S_FF,eps);CHKERRQ(ierr);
485             ierr = MatConvert(S_FF,MATAIJ,MAT_REUSE_MATRIX,&S_FF);CHKERRQ(ierr);
486           }
487           ierr = MatCreateSchurComplement(S_FF,S_FF,S_FE,S_EF,sub_schurs->S_Ej[i],&Stilda_impl);CHKERRQ(ierr);
488           ierr = MatDestroy(&S_FF);CHKERRQ(ierr);
489           ierr = MatDestroy(&S_FE);CHKERRQ(ierr);
490           ierr = MatDestroy(&S_EF);CHKERRQ(ierr);
491           ierr = MatSchurComplementGetKSP(Stilda_impl,&ksp);CHKERRQ(ierr);
492           ierr = KSPSetType(ksp,KSPPREONLY);CHKERRQ(ierr);
493           ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
494           if (is_symmetric) {
495             ierr = PCSetType(pc,PCCHOLESKY);CHKERRQ(ierr);
496           } else {
497             ierr = PCSetType(pc,PCLU);CHKERRQ(ierr);
498           }
499           if (!chop) {
500             ierr = PCFactorSetUseInPlace(pc,PETSC_TRUE);CHKERRQ(ierr);
501           } else {
502             ierr = PCFactorSetMatSolverPackage(pc,MATSOLVERMUMPS);CHKERRQ(ierr);
503           }
504           ierr = KSPSetUp(ksp);CHKERRQ(ierr);
505           ierr = PCBDDCComputeExplicitSchur(Stilda_impl,&Stilda);CHKERRQ(ierr);
506           ierr = MatDestroy(&Stilda_impl);CHKERRQ(ierr);
507         }
508 /*
509         PetscViewerPushFormat(PETSC_VIEWER_STDOUT_SELF,PETSC_VIEWER_ASCII_MATLAB);
510         ierr = MatView(Stilda,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
511         PetscViewerPopFormat(PETSC_VIEWER_STDOUT_SELF);
512 */
513         ierr = MatDestroy(&sub_schurs->S_Ej_tilda[i]);CHKERRQ(ierr);
514         sub_schurs->S_Ej_tilda[i] = Stilda;
515 
516 
517         ierr = MatDenseGetArray(sub_schurs->S_Ej_tilda[i],&vals);CHKERRQ(ierr);
518         if (deluxe) { /* when using deluxe scaling, we need (S_1^-1+S_2^-1)^-1 */
519           PetscScalar *vals2;
520 
521           ierr = MatDenseGetArray(sub_schurs->S_Ej[i],&vals2);CHKERRQ(ierr);
522           /* need to be optimized (cholesky) */
523           ierr = PetscBLASIntCast(subset_size,&B_N);CHKERRQ(ierr);
524           ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
525           PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&B_N,&B_N,vals,&B_N,pivots,&B_ierr));
526           if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRF Lapack routine %d",(int)B_ierr);
527           PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&B_N,vals,&B_N,pivots,work,&B_lwork,&B_ierr));
528           if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRI Lapack routine %d",(int)B_ierr);
529           PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&B_N,&B_N,vals2,&B_N,pivots,&B_ierr));
530           if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRF Lapack routine %d",(int)B_ierr);
531           PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&B_N,vals2,&B_N,pivots,work,&B_lwork,&B_ierr));
532           if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRI Lapack routine %d",(int)B_ierr);
533           ierr = PetscFPTrapPop();CHKERRQ(ierr);
534           ierr = MatSetValues(S_Ej_inv_all,subset_size,dummy_idx,subset_size,dummy_idx,vals2,INSERT_VALUES);CHKERRQ(ierr);
535           ierr = MatDenseRestoreArray(sub_schurs->S_Ej[i],&vals2);CHKERRQ(ierr);
536         }
537         ierr = MatSetValues(S_Ej_tilda_all,subset_size,dummy_idx,subset_size,dummy_idx,vals,INSERT_VALUES);CHKERRQ(ierr);
538         ierr = MatDenseRestoreArray(sub_schurs->S_Ej_tilda[i],&vals);CHKERRQ(ierr);
539       } else {
540         ierr = MatDestroy(&sub_schurs->S_Ej_tilda[i]);CHKERRQ(ierr);
541       }
542       ierr = ISDestroy(&is_E);CHKERRQ(ierr);
543       local_size += subset_size;
544     }
545     ierr = PetscFree(dummy_idx);CHKERRQ(ierr);
546   }
547   ierr = PetscFree(nnz);CHKERRQ(ierr);
548   ierr = MatDestroy(&S_all);CHKERRQ(ierr);
549   ierr = MatAssemblyBegin(sub_schurs->S_Ej_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
550   ierr = MatAssemblyEnd(sub_schurs->S_Ej_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
551   if (compute_Stilda) {
552     ierr = MatAssemblyBegin(S_Ej_tilda_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
553     ierr = MatAssemblyEnd(S_Ej_tilda_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
554     if (deluxe) {
555       ierr = MatAssemblyBegin(S_Ej_inv_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
556       ierr = MatAssemblyEnd(S_Ej_inv_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
557     }
558   }
559 
560   /* Global matrix of all assembled Schur on subsets */
561   ierr = PCBDDCSubsetNumbering(PetscObjectComm((PetscObject)sub_schurs->l2gmap),sub_schurs->l2gmap,local_size,all_local_idx_N+extra,PETSC_NULL,&global_size,&all_local_idx_G);CHKERRQ(ierr);
562   ierr = ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)sub_schurs->l2gmap),1,local_size,all_local_idx_G,PETSC_COPY_VALUES,&l2gmap_subsets);CHKERRQ(ierr);
563   ierr = MatCreateIS(PetscObjectComm((PetscObject)sub_schurs->l2gmap),1,PETSC_DECIDE,PETSC_DECIDE,global_size,global_size,l2gmap_subsets,&work_mat);CHKERRQ(ierr);
564   ierr = ISLocalToGlobalMappingDestroy(&l2gmap_subsets);CHKERRQ(ierr);
565   ierr = MatISSetLocalMat(work_mat,sub_schurs->S_Ej_all);CHKERRQ(ierr);
566   ierr = MatISGetMPIXAIJ(work_mat,MAT_INITIAL_MATRIX,&global_schur_subsets);CHKERRQ(ierr);
567   /* Get local part of (\sum_j S_Ej) */
568   ierr = PetscFree(all_local_idx_N);CHKERRQ(ierr);
569   ierr = ISCreateGeneral(PetscObjectComm((PetscObject)sub_schurs->l2gmap),local_size,all_local_idx_G,PETSC_OWN_POINTER,&temp_is);CHKERRQ(ierr);
570   if (sub_schurs->sum_S_Ej_all) {
571     ierr = MatGetSubMatrixUnsorted(global_schur_subsets,temp_is,temp_is,MAT_REUSE_MATRIX,&sub_schurs->sum_S_Ej_all);CHKERRQ(ierr);
572   } else {
573     ierr = MatGetSubMatrixUnsorted(global_schur_subsets,temp_is,temp_is,MAT_INITIAL_MATRIX,&sub_schurs->sum_S_Ej_all);CHKERRQ(ierr);
574   }
575 
576   if (Stilda_computed) {
577     ierr = MatISSetLocalMat(work_mat,S_Ej_tilda_all);CHKERRQ(ierr);
578     ierr = MatISGetMPIXAIJ(work_mat,MAT_REUSE_MATRIX,&global_schur_subsets);CHKERRQ(ierr);
579     if (sub_schurs->sum_S_Ej_tilda_all) {
580       ierr = MatGetSubMatrixUnsorted(global_schur_subsets,temp_is,temp_is,MAT_REUSE_MATRIX,&sub_schurs->sum_S_Ej_tilda_all);CHKERRQ(ierr);
581     } else {
582       ierr = MatGetSubMatrixUnsorted(global_schur_subsets,temp_is,temp_is,MAT_INITIAL_MATRIX,&sub_schurs->sum_S_Ej_tilda_all);CHKERRQ(ierr);
583     }
584     if (deluxe) {
585       PetscInt    cum;
586       PetscScalar *array,*array2;
587       ierr = MatISSetLocalMat(work_mat,S_Ej_inv_all);CHKERRQ(ierr);
588       ierr = MatISGetMPIXAIJ(work_mat,MAT_REUSE_MATRIX,&global_schur_subsets);CHKERRQ(ierr);
589       if (sub_schurs->sum_S_Ej_inv_all) {
590         ierr = MatGetSubMatrixUnsorted(global_schur_subsets,temp_is,temp_is,MAT_REUSE_MATRIX,&sub_schurs->sum_S_Ej_inv_all);CHKERRQ(ierr);
591       } else {
592         ierr = MatGetSubMatrixUnsorted(global_schur_subsets,temp_is,temp_is,MAT_INITIAL_MATRIX,&sub_schurs->sum_S_Ej_inv_all);CHKERRQ(ierr);
593       }
594       /* invert blocks of sum_S_Ej_inv_all */
595       ierr = MatSeqAIJGetArray(sub_schurs->sum_S_Ej_inv_all,&array);CHKERRQ(ierr);
596       ierr = MatSeqAIJGetArray(sub_schurs->sum_S_Ej_tilda_all,&array2);CHKERRQ(ierr);
597       cum = 0;
598       for (i=0;i<sub_schurs->n_subs;i++) {
599         ierr = ISGetLocalSize(sub_schurs->is_subs[i],&subset_size);CHKERRQ(ierr);
600         if (PetscBTLookup(sub_schurs->computed_Stilda_subs,i)) {
601           /* need to be optimized (cholesky) */
602           ierr = PetscBLASIntCast(subset_size,&B_N);CHKERRQ(ierr);
603           ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
604           PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&B_N,&B_N,array+cum,&B_N,pivots,&B_ierr));
605           if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRF Lapack routine %d",(int)B_ierr);
606           PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&B_N,array+cum,&B_N,pivots,work,&B_lwork,&B_ierr));
607           if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRI Lapack routine %d",(int)B_ierr);
608           PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&B_N,&B_N,array2+cum,&B_N,pivots,&B_ierr));
609           if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRF Lapack routine %d",(int)B_ierr);
610           PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&B_N,array2+cum,&B_N,pivots,work,&B_lwork,&B_ierr));
611           if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRI Lapack routine %d",(int)B_ierr);
612           ierr = PetscFPTrapPop();CHKERRQ(ierr);
613         }
614         cum += subset_size*subset_size;
615       }
616       ierr = MatSeqAIJRestoreArray(sub_schurs->sum_S_Ej_inv_all,&array);CHKERRQ(ierr);
617       ierr = MatSeqAIJRestoreArray(sub_schurs->sum_S_Ej_tilda_all,&array2);CHKERRQ(ierr);
618     }
619   }
620 
621   ierr = PetscFree2(work,pivots);CHKERRQ(ierr);
622   ierr = MatDestroy(&global_schur_subsets);CHKERRQ(ierr);
623   ierr = MatDestroy(&S_Ej_tilda_all);CHKERRQ(ierr);
624   ierr = MatDestroy(&S_Ej_inv_all);CHKERRQ(ierr);
625   ierr = MatDestroy(&work_mat);CHKERRQ(ierr);
626   ierr = ISDestroy(&temp_is);CHKERRQ(ierr);
627   PetscFunctionReturn(0);
628 }
629 
630 #undef __FUNCT__
631 #define __FUNCT__ "PCBDDCSubSchursInit"
632 PetscErrorCode PCBDDCSubSchursInit(PCBDDCSubSchurs sub_schurs, Mat A, Mat S, IS is_I, IS is_B, PCBDDCGraph graph, ISLocalToGlobalMapping BtoNmap, PetscInt seqthreshold)
633 {
634   IS                  *faces,*edges,*all_cc,vertices;
635   PetscInt            *index_sequential,*index_parallel;
636   PetscInt            *auxlocal_sequential,*auxlocal_parallel;
637   PetscInt            *auxglobal_sequential,*auxglobal_parallel;
638   PetscInt            *auxmapping;
639   PetscInt            i,max_subset_size;
640   PetscInt            n_sequential_problems,n_local_sequential_problems,n_parallel_problems,n_local_parallel_problems;
641   PetscInt            n_faces,n_edges,n_all_cc;
642   PetscBool           is_sorted;
643   PetscErrorCode  ierr;
644 
645   PetscFunctionBegin;
646   ierr = ISSorted(is_I,&is_sorted);CHKERRQ(ierr);
647   if (!is_sorted) {
648     SETERRQ(PetscObjectComm((PetscObject)is_I),PETSC_ERR_PLIB,"IS for I dofs should be shorted");
649   }
650   ierr = ISSorted(is_B,&is_sorted);CHKERRQ(ierr);
651   if (!is_sorted) {
652     SETERRQ(PetscObjectComm((PetscObject)is_B),PETSC_ERR_PLIB,"IS for B dofs should be shorted");
653   }
654 
655   /* reset any previous data */
656   ierr = PCBDDCSubSchursReset(sub_schurs);CHKERRQ(ierr);
657 
658   /* get index sets for faces and edges */
659   ierr = PCBDDCGraphGetCandidatesIS(graph,&n_faces,&faces,&n_edges,&edges,&vertices);CHKERRQ(ierr);
660   n_all_cc = n_faces+n_edges;
661   ierr = PetscBTCreate(n_all_cc,&sub_schurs->is_edge);CHKERRQ(ierr);
662   ierr = PetscBTCreate(n_all_cc,&sub_schurs->computed_Stilda_subs);CHKERRQ(ierr);
663   ierr = PetscMalloc1(n_all_cc,&all_cc);CHKERRQ(ierr);
664   for (i=0;i<n_faces;i++) {
665     all_cc[i] = faces[i];
666   }
667   for (i=0;i<n_edges;i++) {
668     all_cc[n_faces+i] = edges[i];
669     ierr = PetscBTSet(sub_schurs->is_edge,n_faces+i);CHKERRQ(ierr);
670   }
671   ierr = PetscFree(faces);CHKERRQ(ierr);
672   ierr = PetscFree(edges);CHKERRQ(ierr);
673 
674   /* map interface's subsets */
675   max_subset_size = 0;
676   for (i=0;i<n_all_cc;i++) {
677     PetscInt subset_size;
678     ierr = ISGetLocalSize(all_cc[i],&subset_size);CHKERRQ(ierr);
679     max_subset_size = PetscMax(max_subset_size,subset_size);
680   }
681   ierr = PetscMalloc1(max_subset_size,&auxmapping);CHKERRQ(ierr);
682   ierr = PetscMalloc2(graph->ncc,&auxlocal_sequential,
683                       graph->ncc,&auxlocal_parallel);CHKERRQ(ierr);
684   ierr = PetscMalloc2(graph->ncc,&index_sequential,
685                       graph->ncc,&index_parallel);CHKERRQ(ierr);
686 
687   /* if threshold is negative uses all sequential problems (possibly using MUMPS) */
688   sub_schurs->use_mumps = PETSC_FALSE;
689   if (seqthreshold < 0) {
690     seqthreshold = max_subset_size;
691 #if defined(PETSC_HAVE_MUMPS)
692     sub_schurs->use_mumps = !!A;
693 #endif
694   }
695 
696 
697   /* determine which problem has to be solved in parallel or sequentially */
698   n_local_sequential_problems = 0;
699   n_local_parallel_problems = 0;
700   for (i=0;i<n_all_cc;i++) {
701     PetscInt       subset_size,j,min_loc = 0;
702     const PetscInt *idxs;
703 
704     ierr = ISGetLocalSize(all_cc[i],&subset_size);CHKERRQ(ierr);
705     ierr = ISGetIndices(all_cc[i],&idxs);CHKERRQ(ierr);
706     ierr = ISLocalToGlobalMappingApply(graph->l2gmap,subset_size,idxs,auxmapping);CHKERRQ(ierr);
707     for (j=1;j<subset_size;j++) {
708       if (auxmapping[j]<auxmapping[min_loc]) {
709         min_loc = j;
710       }
711     }
712     if (subset_size > seqthreshold) {
713       index_parallel[n_local_parallel_problems] = i;
714       auxlocal_parallel[n_local_parallel_problems] = idxs[min_loc];
715       n_local_parallel_problems++;
716     } else {
717       index_sequential[n_local_sequential_problems] = i;
718       auxlocal_sequential[n_local_sequential_problems] = idxs[min_loc];
719       n_local_sequential_problems++;
720     }
721     ierr = ISRestoreIndices(all_cc[i],&idxs);CHKERRQ(ierr);
722   }
723 
724   /* Number parallel problems */
725   auxglobal_parallel = 0;
726   ierr = PCBDDCSubsetNumbering(PetscObjectComm((PetscObject)graph->l2gmap),graph->l2gmap,n_local_parallel_problems,auxlocal_parallel,PETSC_NULL,&n_parallel_problems,&auxglobal_parallel);CHKERRQ(ierr);
727 
728   /* Number sequential problems */
729   auxglobal_sequential = 0;
730   ierr = PCBDDCSubsetNumbering(PetscObjectComm((PetscObject)graph->l2gmap),graph->l2gmap,n_local_sequential_problems,auxlocal_sequential,PETSC_NULL,&n_sequential_problems,&auxglobal_sequential);CHKERRQ(ierr);
731 
732   /* update info in sub_schurs */
733   if (sub_schurs->use_mumps && A) {
734     ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr);
735     sub_schurs->A = A;
736   }
737   ierr = PetscObjectReference((PetscObject)S);CHKERRQ(ierr);
738   sub_schurs->S = S;
739   ierr = PetscObjectReference((PetscObject)is_I);CHKERRQ(ierr);
740   sub_schurs->is_I = is_I;
741   ierr = PetscObjectReference((PetscObject)is_B);CHKERRQ(ierr);
742   sub_schurs->is_B = is_B;
743   ierr = PetscObjectReference((PetscObject)graph->l2gmap);CHKERRQ(ierr);
744   sub_schurs->l2gmap = graph->l2gmap;
745   ierr = PetscObjectReference((PetscObject)BtoNmap);CHKERRQ(ierr);
746   sub_schurs->BtoNmap = BtoNmap;
747   sub_schurs->n_subs_seq = n_local_sequential_problems;
748   sub_schurs->n_subs_par = n_local_parallel_problems;
749   sub_schurs->n_subs_seq_g = n_sequential_problems;
750   sub_schurs->n_subs_par_g = n_parallel_problems;
751   sub_schurs->n_subs = sub_schurs->n_subs_seq + sub_schurs->n_subs_par;
752   sub_schurs->is_subs = all_cc;
753   if (!sub_schurs->use_mumps) { /* for adaptive selection */
754     for (i=0;i<sub_schurs->n_subs;i++) {
755       ierr = ISSort(sub_schurs->is_subs[i]);CHKERRQ(ierr);
756     }
757   }
758   sub_schurs->is_Ej_com = vertices;
759   sub_schurs->index_sequential = index_sequential;
760   sub_schurs->index_parallel = index_parallel;
761   sub_schurs->auxglobal_sequential = auxglobal_sequential;
762   sub_schurs->auxglobal_parallel = auxglobal_parallel;
763 
764 
765   /* allocate space for schur complements */
766   ierr = PetscCalloc2(sub_schurs->n_subs,&sub_schurs->S_Ej,
767                       sub_schurs->n_subs,&sub_schurs->S_Ej_tilda);CHKERRQ(ierr);
768   sub_schurs->S_Ej_all = NULL;
769   sub_schurs->sum_S_Ej_all = NULL;
770   sub_schurs->sum_S_Ej_inv_all = NULL;
771   sub_schurs->sum_S_Ej_tilda_all = NULL;
772   sub_schurs->is_Ej_all = NULL;
773 
774   /* free workspace */
775   ierr = PetscFree(auxmapping);CHKERRQ(ierr);
776   ierr = PetscFree2(auxlocal_sequential,auxlocal_parallel);CHKERRQ(ierr);
777 
778   PetscFunctionReturn(0);
779 }
780 
781 #undef __FUNCT__
782 #define __FUNCT__ "PCBDDCSubSchursCreate"
783 PetscErrorCode PCBDDCSubSchursCreate(PCBDDCSubSchurs *sub_schurs)
784 {
785   PCBDDCSubSchurs schurs_ctx;
786   PetscErrorCode  ierr;
787 
788   PetscFunctionBegin;
789   ierr = PetscNew(&schurs_ctx);CHKERRQ(ierr);
790   schurs_ctx->n_subs = 0;
791   *sub_schurs = schurs_ctx;
792   PetscFunctionReturn(0);
793 }
794 
795 #undef __FUNCT__
796 #define __FUNCT__ "PCBDDCSubSchursDestroy"
797 PetscErrorCode PCBDDCSubSchursDestroy(PCBDDCSubSchurs *sub_schurs)
798 {
799   PetscErrorCode ierr;
800 
801   PetscFunctionBegin;
802   ierr = PCBDDCSubSchursReset(*sub_schurs);CHKERRQ(ierr);
803   ierr = PetscFree(*sub_schurs);CHKERRQ(ierr);
804   PetscFunctionReturn(0);
805 }
806 
807 #undef __FUNCT__
808 #define __FUNCT__ "PCBDDCSubSchursReset"
809 PetscErrorCode PCBDDCSubSchursReset(PCBDDCSubSchurs sub_schurs)
810 {
811   PetscInt       i;
812   PetscErrorCode ierr;
813 
814   PetscFunctionBegin;
815   ierr = MatDestroy(&sub_schurs->A);CHKERRQ(ierr);
816   ierr = MatDestroy(&sub_schurs->S);CHKERRQ(ierr);
817   ierr = ISDestroy(&sub_schurs->is_I);CHKERRQ(ierr);
818   ierr = ISDestroy(&sub_schurs->is_B);CHKERRQ(ierr);
819   ierr = ISLocalToGlobalMappingDestroy(&sub_schurs->l2gmap);CHKERRQ(ierr);
820   ierr = ISLocalToGlobalMappingDestroy(&sub_schurs->BtoNmap);CHKERRQ(ierr);
821   ierr = MatDestroy(&sub_schurs->S_Ej_all);CHKERRQ(ierr);
822   ierr = MatDestroy(&sub_schurs->sum_S_Ej_all);CHKERRQ(ierr);
823   ierr = MatDestroy(&sub_schurs->sum_S_Ej_inv_all);CHKERRQ(ierr);
824   ierr = MatDestroy(&sub_schurs->sum_S_Ej_tilda_all);CHKERRQ(ierr);
825   ierr = ISDestroy(&sub_schurs->is_Ej_all);CHKERRQ(ierr);
826   ierr = ISDestroy(&sub_schurs->is_Ej_com);CHKERRQ(ierr);
827   ierr = PetscBTDestroy(&sub_schurs->is_edge);CHKERRQ(ierr);
828   ierr = PetscBTDestroy(&sub_schurs->computed_Stilda_subs);CHKERRQ(ierr);
829   for (i=0;i<sub_schurs->n_subs;i++) {
830     ierr = ISDestroy(&sub_schurs->is_subs[i]);CHKERRQ(ierr);
831     ierr = MatDestroy(&sub_schurs->S_Ej[i]);CHKERRQ(ierr);
832     ierr = MatDestroy(&sub_schurs->S_Ej_tilda[i]);CHKERRQ(ierr);
833   }
834   ierr = ISDestroy(&sub_schurs->is_I_layer);CHKERRQ(ierr);
835   if (sub_schurs->n_subs) {
836     ierr = PetscFree(sub_schurs->is_subs);CHKERRQ(ierr);
837     ierr = PetscFree2(sub_schurs->S_Ej,sub_schurs->S_Ej_tilda);CHKERRQ(ierr);
838     ierr = PetscFree2(sub_schurs->index_sequential,sub_schurs->index_parallel);CHKERRQ(ierr);
839     ierr = PetscFree(sub_schurs->auxglobal_sequential);CHKERRQ(ierr);
840     ierr = PetscFree(sub_schurs->auxglobal_parallel);CHKERRQ(ierr);
841   }
842   sub_schurs->n_subs = 0;
843   PetscFunctionReturn(0);
844 }
845 
846 #undef __FUNCT__
847 #define __FUNCT__ "PCBDDCAdjGetNextLayer_Private"
848 PETSC_STATIC_INLINE PetscErrorCode PCBDDCAdjGetNextLayer_Private(PetscInt* queue_tip,PetscInt n_prev,PetscBT touched,PetscInt* xadj,PetscInt* adjncy,PetscInt* n_added)
849 {
850   PetscInt       i,j,n;
851   PetscErrorCode ierr;
852 
853   PetscFunctionBegin;
854   n = 0;
855   for (i=-n_prev;i<0;i++) {
856     PetscInt start_dof = queue_tip[i];
857     for (j=xadj[start_dof];j<xadj[start_dof+1];j++) {
858       PetscInt dof = adjncy[j];
859       if (!PetscBTLookup(touched,dof)) {
860         ierr = PetscBTSet(touched,dof);CHKERRQ(ierr);
861         queue_tip[n] = dof;
862         n++;
863       }
864     }
865   }
866   *n_added = n;
867   PetscFunctionReturn(0);
868 }
869