xref: /petsc/src/ksp/pc/impls/bddc/bddcschurs.c (revision bbd98b0ddc3b728de2842e56499debb5c79c81f3) 
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,PetscBool,MatReuse,Mat*);
7 static PetscErrorCode PCBDDCMumpsInteriorSolve(PC,Vec,Vec);
8 static PetscErrorCode PCBDDCMumpsCorrectionSolve(PC,Vec,Vec);
9 
10 #undef __FUNCT__
11 #define __FUNCT__ "PCBDDCMumpsCorrectionSolve_Private"
12 static PetscErrorCode PCBDDCMumpsCorrectionSolve_Private(PC pc, Vec rhs, Vec sol, PetscBool transpose)
13 {
14   PCBDDCReuseMumps ctx;
15 #if defined(PETSC_HAVE_MUMPS)
16   PetscInt         ival;
17 #endif
18   PetscErrorCode   ierr;
19 
20   PetscFunctionBegin;
21   ierr = PCShellGetContext(pc,(void **)&ctx);CHKERRQ(ierr);
22 #if defined(PETSC_HAVE_MUMPS)
23   ierr = MatMumpsGetIcntl(ctx->F,26,&ival);CHKERRQ(ierr);
24   ierr = MatMumpsSetIcntl(ctx->F,26,-1);CHKERRQ(ierr);
25 #endif
26   if (transpose) {
27     ierr = MatSolveTranspose(ctx->F,rhs,sol);CHKERRQ(ierr);
28   } else {
29     ierr = MatSolve(ctx->F,rhs,sol);CHKERRQ(ierr);
30   }
31 #if defined(PETSC_HAVE_MUMPS)
32   ierr = MatMumpsSetIcntl(ctx->F,26,ival);CHKERRQ(ierr);
33 #endif
34   PetscFunctionReturn(0);
35 }
36 
37 #undef __FUNCT__
38 #define __FUNCT__ "PCBDDCMumpsCorrectionSolve"
39 static PetscErrorCode PCBDDCMumpsCorrectionSolve(PC pc, Vec rhs, Vec sol)
40 {
41   PetscErrorCode   ierr;
42 
43   PetscFunctionBegin;
44   ierr = PCBDDCMumpsCorrectionSolve_Private(pc,rhs,sol,PETSC_FALSE);CHKERRQ(ierr);
45   PetscFunctionReturn(0);
46 }
47 
48 #undef __FUNCT__
49 #define __FUNCT__ "PCBDDCMumpsCorrectionSolveTranspose"
50 static PetscErrorCode PCBDDCMumpsCorrectionSolveTranspose(PC pc, Vec rhs, Vec sol)
51 {
52   PetscErrorCode   ierr;
53 
54   PetscFunctionBegin;
55   ierr = PCBDDCMumpsCorrectionSolve_Private(pc,rhs,sol,PETSC_TRUE);CHKERRQ(ierr);
56   PetscFunctionReturn(0);
57 }
58 
59 #undef __FUNCT__
60 #define __FUNCT__ "PCBDDCReuseMumpsReset"
61 static PetscErrorCode PCBDDCReuseMumpsReset(PCBDDCReuseMumps reuse)
62 {
63   PetscErrorCode ierr;
64 
65   PetscFunctionBegin;
66   ierr = MatDestroy(&reuse->F);CHKERRQ(ierr);
67   ierr = VecDestroy(&reuse->sol);CHKERRQ(ierr);
68   ierr = VecDestroy(&reuse->rhs);CHKERRQ(ierr);
69   ierr = PCDestroy(&reuse->interior_solver);CHKERRQ(ierr);
70   ierr = PCDestroy(&reuse->correction_solver);CHKERRQ(ierr);
71   ierr = ISDestroy(&reuse->is_R);CHKERRQ(ierr);
72   ierr = ISDestroy(&reuse->is_B);CHKERRQ(ierr);
73   ierr = VecScatterDestroy(&reuse->correction_scatter_B);CHKERRQ(ierr);
74   ierr = VecDestroy(&reuse->sol_B);CHKERRQ(ierr);
75   ierr = VecDestroy(&reuse->rhs_B);CHKERRQ(ierr);
76   PetscFunctionReturn(0);
77 }
78 
79 #undef __FUNCT__
80 #define __FUNCT__ "PCBDDCMumpsInteriorSolve_Private"
81 static PetscErrorCode PCBDDCMumpsInteriorSolve_Private(PC pc, Vec rhs, Vec sol, PetscBool transpose)
82 {
83   PCBDDCReuseMumps ctx;
84   PetscScalar      *array,*array_mumps;
85 #if defined(PETSC_HAVE_MUMPS)
86   PetscInt         ival;
87 #endif
88   PetscErrorCode   ierr;
89 
90   PetscFunctionBegin;
91   ierr = PCShellGetContext(pc,(void **)&ctx);CHKERRQ(ierr);
92 #if defined(PETSC_HAVE_MUMPS)
93   ierr = MatMumpsGetIcntl(ctx->F,26,&ival);CHKERRQ(ierr);
94   ierr = MatMumpsSetIcntl(ctx->F,26,0);CHKERRQ(ierr);
95 #endif
96   /* copy rhs into factored matrix workspace (can it be avoided?, MatSolve_MUMPS has another copy b->x internally) */
97   ierr = VecGetArrayRead(rhs,(const PetscScalar**)&array);CHKERRQ(ierr);
98   ierr = VecGetArray(ctx->rhs,&array_mumps);CHKERRQ(ierr);
99   ierr = PetscMemcpy(array_mumps,array,ctx->n_I*sizeof(PetscScalar));CHKERRQ(ierr);
100   ierr = VecRestoreArray(ctx->rhs,&array_mumps);CHKERRQ(ierr);
101   ierr = VecRestoreArrayRead(rhs,(const PetscScalar**)&array);CHKERRQ(ierr);
102 
103   if (transpose) {
104     ierr = MatSolveTranspose(ctx->F,ctx->rhs,ctx->sol);CHKERRQ(ierr);
105   } else {
106     ierr = MatSolve(ctx->F,ctx->rhs,ctx->sol);CHKERRQ(ierr);
107   }
108 
109   /* get back data to caller worskpace */
110   ierr = VecGetArrayRead(ctx->sol,(const PetscScalar**)&array_mumps);CHKERRQ(ierr);
111   ierr = VecGetArray(sol,&array);CHKERRQ(ierr);
112   ierr = PetscMemcpy(array,array_mumps,ctx->n_I*sizeof(PetscScalar));CHKERRQ(ierr);
113   ierr = VecRestoreArray(sol,&array);CHKERRQ(ierr);
114   ierr = VecRestoreArrayRead(ctx->sol,(const PetscScalar**)&array_mumps);CHKERRQ(ierr);
115 #if defined(PETSC_HAVE_MUMPS)
116   ierr = MatMumpsSetIcntl(ctx->F,26,ival);CHKERRQ(ierr);
117 #endif
118   PetscFunctionReturn(0);
119 }
120 
121 #undef __FUNCT__
122 #define __FUNCT__ "PCBDDCMumpsInteriorSolve"
123 static PetscErrorCode PCBDDCMumpsInteriorSolve(PC pc, Vec rhs, Vec sol)
124 {
125   PetscErrorCode   ierr;
126 
127   PetscFunctionBegin;
128   ierr = PCBDDCMumpsInteriorSolve_Private(pc,rhs,sol,PETSC_FALSE);CHKERRQ(ierr);
129   PetscFunctionReturn(0);
130 }
131 
132 #undef __FUNCT__
133 #define __FUNCT__ "PCBDDCMumpsInteriorSolveTranspose"
134 static PetscErrorCode PCBDDCMumpsInteriorSolveTranspose(PC pc, Vec rhs, Vec sol)
135 {
136   PetscErrorCode   ierr;
137 
138   PetscFunctionBegin;
139   ierr = PCBDDCMumpsInteriorSolve_Private(pc,rhs,sol,PETSC_TRUE);CHKERRQ(ierr);
140   PetscFunctionReturn(0);
141 }
142 
143 #undef __FUNCT__
144 #define __FUNCT__ "PCBDDCComputeExplicitSchur"
145 static PetscErrorCode PCBDDCComputeExplicitSchur(Mat M, PetscBool issym, MatReuse reuse, Mat *S)
146 {
147   Mat            B, C, D, Bd, Cd, AinvBd;
148   KSP            ksp;
149   PC             pc;
150   PetscBool      isLU, isILU, isCHOL, Bdense, Cdense;
151   PetscReal      fill = 2.0;
152   PetscInt       n_I;
153   PetscMPIInt    size;
154   PetscErrorCode ierr;
155 
156   PetscFunctionBegin;
157   ierr = MPI_Comm_size(PetscObjectComm((PetscObject)M),&size);CHKERRQ(ierr);
158   if (size != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for parallel matrices");
159   if (reuse == MAT_REUSE_MATRIX) {
160     PetscBool Sdense;
161 
162     ierr = PetscObjectTypeCompare((PetscObject)*S, MATSEQDENSE, &Sdense);CHKERRQ(ierr);
163     if (!Sdense) SETERRQ(PetscObjectComm((PetscObject)M),PETSC_ERR_SUP,"S should dense");
164   }
165   ierr = MatSchurComplementGetSubMatrices(M, NULL, NULL, &B, &C, &D);CHKERRQ(ierr);
166   ierr = MatSchurComplementGetKSP(M, &ksp);CHKERRQ(ierr);
167   ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr);
168   ierr = PetscObjectTypeCompare((PetscObject) pc, PCLU, &isLU);CHKERRQ(ierr);
169   ierr = PetscObjectTypeCompare((PetscObject) pc, PCILU, &isILU);CHKERRQ(ierr);
170   ierr = PetscObjectTypeCompare((PetscObject) pc, PCCHOLESKY, &isCHOL);CHKERRQ(ierr);
171   ierr = PetscObjectTypeCompare((PetscObject) B, MATSEQDENSE, &Bdense);CHKERRQ(ierr);
172   ierr = PetscObjectTypeCompare((PetscObject) C, MATSEQDENSE, &Cdense);CHKERRQ(ierr);
173   ierr = MatGetSize(B,&n_I,NULL);CHKERRQ(ierr);
174   if (n_I) {
175     if (!Bdense) {
176       ierr = MatConvert(B, MATSEQDENSE, MAT_INITIAL_MATRIX, &Bd);CHKERRQ(ierr);
177     } else {
178       Bd = B;
179     }
180 
181     if (isLU || isILU || isCHOL) {
182       Mat fact;
183       ierr = KSPSetUp(ksp);CHKERRQ(ierr);
184       ierr = PCFactorGetMatrix(pc, &fact);CHKERRQ(ierr);
185       ierr = MatDuplicate(Bd, MAT_DO_NOT_COPY_VALUES, &AinvBd);CHKERRQ(ierr);
186       ierr = MatMatSolve(fact, Bd, AinvBd);CHKERRQ(ierr);
187     } else {
188       PetscBool ex = PETSC_TRUE;
189 
190       if (ex) {
191         Mat Ainvd;
192 
193         ierr = PCComputeExplicitOperator(pc, &Ainvd);CHKERRQ(ierr);
194         ierr = MatMatMult(Ainvd, Bd, MAT_INITIAL_MATRIX, fill, &AinvBd);CHKERRQ(ierr);
195         ierr = MatDestroy(&Ainvd);CHKERRQ(ierr);
196       } else {
197         Vec         sol,rhs;
198         PetscScalar *arrayrhs,*arraysol;
199         PetscInt    i,nrhs,n;
200 
201         ierr = MatDuplicate(Bd, MAT_DO_NOT_COPY_VALUES, &AinvBd);CHKERRQ(ierr);
202         ierr = MatGetSize(Bd,&n,&nrhs);CHKERRQ(ierr);
203         ierr = MatDenseGetArray(Bd,&arrayrhs);CHKERRQ(ierr);
204         ierr = MatDenseGetArray(AinvBd,&arraysol);CHKERRQ(ierr);
205         ierr = KSPGetSolution(ksp,&sol);CHKERRQ(ierr);
206         ierr = KSPGetRhs(ksp,&rhs);CHKERRQ(ierr);
207         for (i=0;i<nrhs;i++) {
208           ierr = VecPlaceArray(rhs,arrayrhs+i*n);CHKERRQ(ierr);
209           ierr = VecPlaceArray(sol,arraysol+i*n);CHKERRQ(ierr);
210           ierr = KSPSolve(ksp,rhs,sol);CHKERRQ(ierr);
211           ierr = VecResetArray(rhs);CHKERRQ(ierr);
212           ierr = VecResetArray(sol);CHKERRQ(ierr);
213         }
214         ierr = MatDenseRestoreArray(Bd,&arrayrhs);CHKERRQ(ierr);
215         ierr = MatDenseRestoreArray(AinvBd,&arrayrhs);CHKERRQ(ierr);
216       }
217     }
218     if (!Bdense & !issym) {
219       ierr = MatDestroy(&Bd);CHKERRQ(ierr);
220     }
221 
222     if (!issym) {
223       if (!Cdense) {
224         ierr = MatConvert(C, MATSEQDENSE, MAT_INITIAL_MATRIX, &Cd);CHKERRQ(ierr);
225       } else {
226         Cd = C;
227       }
228       ierr = MatMatMult(Cd, AinvBd, reuse, fill, S);CHKERRQ(ierr);
229       if (!Cdense) {
230         ierr = MatDestroy(&Cd);CHKERRQ(ierr);
231       }
232     } else {
233       ierr = MatTransposeMatMult(Bd, AinvBd, reuse, fill, S);CHKERRQ(ierr);
234       if (!Bdense) {
235         ierr = MatDestroy(&Bd);CHKERRQ(ierr);
236       }
237     }
238     ierr = MatDestroy(&AinvBd);CHKERRQ(ierr);
239   }
240 
241   if (D) {
242     Mat       Dd;
243     PetscBool Ddense;
244 
245     ierr = PetscObjectTypeCompare((PetscObject)D,MATSEQDENSE,&Ddense);CHKERRQ(ierr);
246     if (!Ddense) {
247       ierr = MatConvert(D, MATSEQDENSE, MAT_INITIAL_MATRIX, &Dd);CHKERRQ(ierr);
248     } else {
249       Dd = D;
250     }
251     if (n_I) {
252       ierr = MatAYPX(*S,-1.0,Dd,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
253     } else {
254       if (reuse == MAT_INITIAL_MATRIX) {
255         ierr = MatDuplicate(Dd,MAT_COPY_VALUES,S);CHKERRQ(ierr);
256       } else {
257         ierr = MatCopy(Dd,*S,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
258       }
259     }
260     if (!Ddense) {
261       ierr = MatDestroy(&Dd);CHKERRQ(ierr);
262     }
263   } else {
264     ierr = MatScale(*S,-1.0);CHKERRQ(ierr);
265   }
266   PetscFunctionReturn(0);
267 }
268 
269 #undef __FUNCT__
270 #define __FUNCT__ "PCBDDCSubSchursSetUp"
271 PetscErrorCode PCBDDCSubSchursSetUp(PCBDDCSubSchurs sub_schurs, Mat Ain, Mat Sin, PetscInt xadj[], PetscInt adjncy[], PetscInt nlayers, PetscBool faster_deluxe, PetscBool compute_Stilda, PetscBool reuse_solvers)
272 {
273   Mat                    F,A_II,A_IB,A_BI,A_BB,AE_II;
274   Mat                    S_all;
275   Mat                    global_schur_subsets,work_mat,*submats;
276   ISLocalToGlobalMapping l2gmap_subsets;
277   IS                     is_I,is_I_layer;
278   IS                     all_subsets,all_subsets_mult,all_subsets_n;
279   PetscInt               *nnz,*all_local_idx_N;
280   PetscInt               *auxnum1,*auxnum2;
281   PetscInt               i,subset_size,max_subset_size;
282   PetscInt               extra,local_size,global_size;
283   PetscBLASInt           B_N,B_ierr,B_lwork,*pivots;
284   PetscScalar            *Bwork;
285   PetscSubcomm           subcomm;
286   PetscMPIInt            color,rank;
287   MPI_Comm               comm_n;
288   PetscErrorCode         ierr;
289 
290   PetscFunctionBegin;
291   /* update info in sub_schurs */
292   ierr = MatDestroy(&sub_schurs->A);CHKERRQ(ierr);
293   ierr = MatDestroy(&sub_schurs->S);CHKERRQ(ierr);
294   sub_schurs->is_hermitian = PETSC_FALSE;
295   sub_schurs->is_posdef = PETSC_FALSE;
296   if (Ain) {
297     PetscBool isseqaij;
298     /* determine if we are dealing with hermitian positive definite problems */
299 #if !defined(PETSC_USE_COMPLEX)
300     if (Ain->symmetric_set) {
301       sub_schurs->is_hermitian = Ain->symmetric;
302     }
303 #else
304     if (Ain->hermitian_set) {
305       sub_schurs->is_hermitian = Ain->hermitian;
306     }
307 #endif
308     if (Ain->spd_set) {
309       sub_schurs->is_posdef = Ain->spd;
310     }
311 
312     /* check */
313     ierr = PetscObjectTypeCompare((PetscObject)Ain,MATSEQAIJ,&isseqaij);CHKERRQ(ierr);
314     if (compute_Stilda && (!sub_schurs->is_hermitian || !sub_schurs->is_posdef)) {
315       PetscInt lsize;
316 
317       ierr = MatGetSize(Ain,&lsize,NULL);CHKERRQ(ierr);
318       if (lsize) {
319         PetscScalar val;
320         PetscReal   norm;
321         Vec         vec1,vec2,vec3;
322 
323         ierr = MatCreateVecs(Ain,&vec1,&vec2);CHKERRQ(ierr);
324         ierr = VecDuplicate(vec1,&vec3);CHKERRQ(ierr);
325         ierr = VecSetRandom(vec1,NULL);CHKERRQ(ierr);
326         ierr = MatMult(Ain,vec1,vec2);CHKERRQ(ierr);
327 #if !defined(PETSC_USE_COMPLEX)
328         ierr = MatMultTranspose(Ain,vec1,vec3);CHKERRQ(ierr);
329 #else
330         ierr = MatMultHermitianTranspose(Ain,vec1,vec3);CHKERRQ(ierr);
331 #endif
332         ierr = VecAXPY(vec3,-1.0,vec2);CHKERRQ(ierr);
333         ierr = VecNorm(vec3,NORM_INFINITY,&norm);CHKERRQ(ierr);
334         if (norm > PetscSqrtReal(PETSC_SMALL)) {
335           sub_schurs->is_hermitian = PETSC_FALSE;
336         } else {
337           sub_schurs->is_hermitian = PETSC_TRUE;
338         }
339         ierr = VecDot(vec1,vec2,&val);CHKERRQ(ierr);
340         if (PetscRealPart(val) > 0. && PetscAbsReal(PetscImaginaryPart(val)) < PETSC_SMALL) sub_schurs->is_posdef = PETSC_TRUE;
341         ierr = VecDestroy(&vec1);CHKERRQ(ierr);
342         ierr = VecDestroy(&vec2);CHKERRQ(ierr);
343         ierr = VecDestroy(&vec3);CHKERRQ(ierr);
344       } else {
345         sub_schurs->is_hermitian = PETSC_TRUE;
346         sub_schurs->is_posdef = PETSC_TRUE;
347       }
348     }
349     if (isseqaij) {
350       ierr = PetscObjectReference((PetscObject)Ain);CHKERRQ(ierr);
351       sub_schurs->A = Ain;
352     } else {
353       ierr = MatConvert(Ain,MATSEQAIJ,MAT_INITIAL_MATRIX,&sub_schurs->A);CHKERRQ(ierr);
354     }
355   }
356   if (compute_Stilda && (!sub_schurs->is_hermitian || !sub_schurs->is_posdef)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"General matrix pencils are not currently supported (%D,%D)",sub_schurs->is_hermitian,sub_schurs->is_posdef);
357 
358   ierr = PetscObjectReference((PetscObject)Sin);CHKERRQ(ierr);
359   sub_schurs->S = Sin;
360   if (sub_schurs->use_mumps) {
361     sub_schurs->use_mumps = (PetscBool)(!!sub_schurs->A);
362   }
363 
364   /* preliminary checks */
365   if (!sub_schurs->use_mumps && compute_Stilda) SETERRQ(PetscObjectComm((PetscObject)sub_schurs->l2gmap),PETSC_ERR_SUP,"Adaptive selection of constraints requires MUMPS");
366 
367   /* restrict work on active processes */
368   color = 0;
369   if (!sub_schurs->n_subs) color = 1; /* this can happen if we are in a multilevel case or if the subdomain is disconnected */
370   ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)sub_schurs->l2gmap),&rank);CHKERRQ(ierr);
371   ierr = PetscSubcommCreate(PetscObjectComm((PetscObject)sub_schurs->l2gmap),&subcomm);CHKERRQ(ierr);
372   ierr = PetscSubcommSetNumber(subcomm,2);CHKERRQ(ierr);
373   ierr = PetscSubcommSetTypeGeneral(subcomm,color,rank);CHKERRQ(ierr);
374   ierr = PetscCommDuplicate(PetscSubcommChild(subcomm),&comm_n,NULL);CHKERRQ(ierr);
375   ierr = PetscSubcommDestroy(&subcomm);CHKERRQ(ierr);
376   if (!sub_schurs->n_subs) {
377     ierr = PetscCommDestroy(&comm_n);CHKERRQ(ierr);
378     PetscFunctionReturn(0);
379   }
380   /* ierr = PetscCommDuplicate(PetscObjectComm((PetscObject)sub_schurs->l2gmap),&comm_n,NULL);CHKERRQ(ierr); */
381 
382   /* get Schur complement matrices */
383   if (!sub_schurs->use_mumps) {
384     Mat       tA_IB,tA_BI,tA_BB;
385     PetscBool isseqsbaij;
386     ierr = MatSchurComplementGetSubMatrices(sub_schurs->S,&A_II,NULL,&tA_IB,&tA_BI,&tA_BB);CHKERRQ(ierr);
387     ierr = PetscObjectTypeCompare((PetscObject)tA_BB,MATSEQSBAIJ,&isseqsbaij);CHKERRQ(ierr);
388     if (isseqsbaij) {
389       ierr = MatConvert(tA_BB,MATSEQAIJ,MAT_INITIAL_MATRIX,&A_BB);CHKERRQ(ierr);
390       ierr = MatConvert(tA_IB,MATSEQAIJ,MAT_INITIAL_MATRIX,&A_IB);CHKERRQ(ierr);
391       ierr = MatConvert(tA_BI,MATSEQAIJ,MAT_INITIAL_MATRIX,&A_BI);CHKERRQ(ierr);
392     } else {
393       ierr = PetscObjectReference((PetscObject)tA_BB);CHKERRQ(ierr);
394       A_BB = tA_BB;
395       ierr = PetscObjectReference((PetscObject)tA_IB);CHKERRQ(ierr);
396       A_IB = tA_IB;
397       ierr = PetscObjectReference((PetscObject)tA_BI);CHKERRQ(ierr);
398       A_BI = tA_BI;
399     }
400   } else {
401     A_II = NULL;
402     A_IB = NULL;
403     A_BI = NULL;
404     A_BB = NULL;
405   }
406   S_all = NULL;
407 
408   /* determine interior problems */
409   ierr = ISGetLocalSize(sub_schurs->is_I,&i);CHKERRQ(ierr);
410   if (nlayers >= 0 && i) { /* Interior problems can be different from the original one */
411     PetscBT                touched;
412     const PetscInt*        idx_B;
413     PetscInt               n_I,n_B,n_local_dofs,n_prev_added,j,layer,*local_numbering;
414 
415     if (!xadj || !adjncy) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Cannot request layering without adjacency");
416     /* get sizes */
417     ierr = ISGetLocalSize(sub_schurs->is_I,&n_I);CHKERRQ(ierr);
418     ierr = ISGetLocalSize(sub_schurs->is_B,&n_B);CHKERRQ(ierr);
419 
420     ierr = PetscMalloc1(n_I+n_B,&local_numbering);CHKERRQ(ierr);
421     ierr = PetscBTCreate(n_I+n_B,&touched);CHKERRQ(ierr);
422     ierr = PetscBTMemzero(n_I+n_B,touched);CHKERRQ(ierr);
423 
424     /* all boundary dofs must be skipped when adding layers */
425     ierr = ISGetIndices(sub_schurs->is_B,&idx_B);CHKERRQ(ierr);
426     for (j=0;j<n_B;j++) {
427       ierr = PetscBTSet(touched,idx_B[j]);CHKERRQ(ierr);
428     }
429     ierr = PetscMemcpy(local_numbering,idx_B,n_B*sizeof(PetscInt));CHKERRQ(ierr);
430     ierr = ISRestoreIndices(sub_schurs->is_B,&idx_B);CHKERRQ(ierr);
431 
432     /* add prescribed number of layers of dofs */
433     n_local_dofs = n_B;
434     n_prev_added = n_B;
435     for (layer=0;layer<nlayers;layer++) {
436       PetscInt n_added;
437       if (n_local_dofs == n_I+n_B) break;
438       if (n_local_dofs > n_I+n_B) 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);
439       ierr = PCBDDCAdjGetNextLayer_Private(local_numbering+n_local_dofs,n_prev_added,touched,xadj,adjncy,&n_added);CHKERRQ(ierr);
440       n_prev_added = n_added;
441       n_local_dofs += n_added;
442       if (!n_added) break;
443     }
444     ierr = PetscBTDestroy(&touched);CHKERRQ(ierr);
445 
446     /* IS for I layer dofs in original numbering */
447     ierr = ISCreateGeneral(PetscObjectComm((PetscObject)sub_schurs->is_I),n_local_dofs-n_B,local_numbering+n_B,PETSC_COPY_VALUES,&is_I_layer);CHKERRQ(ierr);
448     ierr = PetscFree(local_numbering);CHKERRQ(ierr);
449     ierr = ISSort(is_I_layer);CHKERRQ(ierr);
450     /* IS for I layer dofs in I numbering */
451     if (!sub_schurs->use_mumps) {
452       ISLocalToGlobalMapping ItoNmap;
453       ierr = ISLocalToGlobalMappingCreateIS(sub_schurs->is_I,&ItoNmap);CHKERRQ(ierr);
454       ierr = ISGlobalToLocalMappingApplyIS(ItoNmap,IS_GTOLM_DROP,is_I_layer,&is_I);CHKERRQ(ierr);
455       ierr = ISLocalToGlobalMappingDestroy(&ItoNmap);CHKERRQ(ierr);
456 
457       /* II block */
458       ierr = MatGetSubMatrix(A_II,is_I,is_I,MAT_INITIAL_MATRIX,&AE_II);CHKERRQ(ierr);
459     }
460   } else {
461     PetscInt n_I;
462 
463     /* IS for I dofs in original numbering */
464     ierr = PetscObjectReference((PetscObject)sub_schurs->is_I);CHKERRQ(ierr);
465     is_I_layer = sub_schurs->is_I;
466 
467     /* IS for I dofs in I numbering (strided 1) */
468     if (!sub_schurs->use_mumps) {
469       ierr = ISGetSize(sub_schurs->is_I,&n_I);CHKERRQ(ierr);
470       ierr = ISCreateStride(PetscObjectComm((PetscObject)sub_schurs->is_I),n_I,0,1,&is_I);CHKERRQ(ierr);
471 
472       /* II block is the same */
473       ierr = PetscObjectReference((PetscObject)A_II);CHKERRQ(ierr);
474       AE_II = A_II;
475     }
476   }
477 
478   /* Get info on subset sizes and sum of all subsets sizes */
479   max_subset_size = 0;
480   local_size = 0;
481   for (i=0;i<sub_schurs->n_subs;i++) {
482     ierr = ISGetLocalSize(sub_schurs->is_subs[i],&subset_size);CHKERRQ(ierr);
483     max_subset_size = PetscMax(subset_size,max_subset_size);
484     local_size += subset_size;
485   }
486 
487   /* Work arrays for local indices */
488   extra = 0;
489   if (sub_schurs->use_mumps && is_I_layer) {
490     ierr = ISGetLocalSize(is_I_layer,&extra);CHKERRQ(ierr);
491   }
492   ierr = PetscMalloc1(local_size+extra,&all_local_idx_N);CHKERRQ(ierr);
493   if (extra) {
494     const PetscInt *idxs;
495     ierr = ISGetIndices(is_I_layer,&idxs);CHKERRQ(ierr);
496     ierr = PetscMemcpy(all_local_idx_N,idxs,extra*sizeof(PetscInt));CHKERRQ(ierr);
497     ierr = ISRestoreIndices(is_I_layer,&idxs);CHKERRQ(ierr);
498   }
499   ierr = PetscMalloc1(local_size,&nnz);CHKERRQ(ierr);
500   ierr = PetscMalloc1(sub_schurs->n_subs,&auxnum1);CHKERRQ(ierr);
501   ierr = PetscMalloc1(sub_schurs->n_subs,&auxnum2);CHKERRQ(ierr);
502 
503   /* Get local indices in local numbering */
504   local_size = 0;
505   for (i=0;i<sub_schurs->n_subs;i++) {
506     PetscInt j;
507     const    PetscInt *idxs;
508 
509     ierr = ISGetLocalSize(sub_schurs->is_subs[i],&subset_size);CHKERRQ(ierr);
510     ierr = ISGetIndices(sub_schurs->is_subs[i],&idxs);CHKERRQ(ierr);
511     /* start (smallest in global ordering) and multiplicity */
512     auxnum1[i] = idxs[0];
513     auxnum2[i] = subset_size;
514     /* subset indices in local numbering */
515     ierr = PetscMemcpy(all_local_idx_N+local_size+extra,idxs,subset_size*sizeof(PetscInt));CHKERRQ(ierr);
516     ierr = ISRestoreIndices(sub_schurs->is_subs[i],&idxs);CHKERRQ(ierr);
517     for (j=0;j<subset_size;j++) nnz[local_size+j] = subset_size;
518     local_size += subset_size;
519   }
520 
521   /* allocate extra workspace needed only for GETRI */
522   Bwork = NULL;
523   pivots = NULL;
524   if (local_size && (!sub_schurs->is_hermitian || !sub_schurs->is_posdef)) {
525     PetscScalar lwork;
526 
527     B_lwork = -1;
528     ierr = PetscBLASIntCast(local_size,&B_N);CHKERRQ(ierr);
529     ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
530     PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&B_N,Bwork,&B_N,pivots,&lwork,&B_lwork,&B_ierr));
531     ierr = PetscFPTrapPop();CHKERRQ(ierr);
532     if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in query to GETRI Lapack routine %d",(int)B_ierr);
533     ierr = PetscBLASIntCast((PetscInt)PetscRealPart(lwork),&B_lwork);CHKERRQ(ierr);
534     ierr = PetscMalloc2(B_lwork,&Bwork,B_N,&pivots);CHKERRQ(ierr);
535   }
536 
537   /* prepare parallel matrices for summing up properly schurs on subsets */
538   ierr = ISCreateGeneral(comm_n,sub_schurs->n_subs,auxnum1,PETSC_OWN_POINTER,&all_subsets_n);CHKERRQ(ierr);
539   ierr = ISLocalToGlobalMappingApplyIS(sub_schurs->l2gmap,all_subsets_n,&all_subsets);CHKERRQ(ierr);
540   ierr = ISDestroy(&all_subsets_n);CHKERRQ(ierr);
541   ierr = ISCreateGeneral(comm_n,sub_schurs->n_subs,auxnum2,PETSC_OWN_POINTER,&all_subsets_mult);CHKERRQ(ierr);
542   ierr = PCBDDCSubsetNumbering(all_subsets,all_subsets_mult,&global_size,&all_subsets_n);CHKERRQ(ierr);
543   ierr = ISDestroy(&all_subsets);CHKERRQ(ierr);
544   ierr = ISDestroy(&all_subsets_mult);CHKERRQ(ierr);
545   ierr = ISGetLocalSize(all_subsets_n,&i);CHKERRQ(ierr);
546   if (i != local_size) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Invalid size of new subset! %D != %D",i,local_size);
547   ierr = ISLocalToGlobalMappingCreateIS(all_subsets_n,&l2gmap_subsets);CHKERRQ(ierr);
548   ierr = MatCreateIS(comm_n,1,PETSC_DECIDE,PETSC_DECIDE,global_size,global_size,l2gmap_subsets,NULL,&work_mat);CHKERRQ(ierr);
549   ierr = ISLocalToGlobalMappingDestroy(&l2gmap_subsets);CHKERRQ(ierr);
550   ierr = MatCreate(PetscObjectComm((PetscObject)work_mat),&global_schur_subsets);CHKERRQ(ierr);
551   ierr = MatSetSizes(global_schur_subsets,PETSC_DECIDE,PETSC_DECIDE,global_size,global_size);CHKERRQ(ierr);
552   ierr = MatSetType(global_schur_subsets,MATMPIAIJ);CHKERRQ(ierr);
553 
554   /* subset indices in local boundary numbering */
555   if (!sub_schurs->is_Ej_all) {
556     PetscInt *all_local_idx_B;
557 
558     ierr = PetscMalloc1(local_size,&all_local_idx_B);CHKERRQ(ierr);
559     ierr = ISGlobalToLocalMappingApply(sub_schurs->BtoNmap,IS_GTOLM_DROP,local_size,all_local_idx_N+extra,&subset_size,all_local_idx_B);CHKERRQ(ierr);
560     if (subset_size != local_size) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in sub_schurs serial (BtoNmap)! %D != %D\n",subset_size,local_size);
561     ierr = ISCreateGeneral(PETSC_COMM_SELF,local_size,all_local_idx_B,PETSC_OWN_POINTER,&sub_schurs->is_Ej_all);CHKERRQ(ierr);
562   }
563 
564   /* Local matrix of all local Schur on subsets (transposed) */
565   if (!sub_schurs->S_Ej_all) {
566     ierr = MatCreate(PETSC_COMM_SELF,&sub_schurs->S_Ej_all);CHKERRQ(ierr);
567     ierr = MatSetSizes(sub_schurs->S_Ej_all,PETSC_DECIDE,PETSC_DECIDE,local_size,local_size);CHKERRQ(ierr);
568     ierr = MatSetType(sub_schurs->S_Ej_all,MATAIJ);CHKERRQ(ierr);
569     ierr = MatSeqAIJSetPreallocation(sub_schurs->S_Ej_all,0,nnz);CHKERRQ(ierr);
570   }
571 
572   /* Compute Schur complements explicitly */
573   F = NULL;
574   if (!sub_schurs->use_mumps) {
575     Mat         S_Ej_expl;
576     PetscScalar *work;
577     PetscInt    j,*dummy_idx;
578     PetscBool   Sdense;
579 
580     ierr = PetscMalloc2(max_subset_size,&dummy_idx,max_subset_size*max_subset_size,&work);CHKERRQ(ierr);
581     local_size = 0;
582     for (i=0;i<sub_schurs->n_subs;i++) {
583       IS  is_subset_B;
584       Mat AE_EE,AE_IE,AE_EI,S_Ej;
585 
586       /* subsets in original and boundary numbering */
587       ierr = ISGlobalToLocalMappingApplyIS(sub_schurs->BtoNmap,IS_GTOLM_DROP,sub_schurs->is_subs[i],&is_subset_B);CHKERRQ(ierr);
588       /* EE block */
589       ierr = MatGetSubMatrix(A_BB,is_subset_B,is_subset_B,MAT_INITIAL_MATRIX,&AE_EE);CHKERRQ(ierr);
590       /* IE block */
591       ierr = MatGetSubMatrix(A_IB,is_I,is_subset_B,MAT_INITIAL_MATRIX,&AE_IE);CHKERRQ(ierr);
592       /* EI block */
593       if (sub_schurs->is_hermitian) {
594         ierr = MatCreateTranspose(AE_IE,&AE_EI);CHKERRQ(ierr);
595       } else {
596         ierr = MatGetSubMatrix(A_BI,is_subset_B,is_I,MAT_INITIAL_MATRIX,&AE_EI);CHKERRQ(ierr);
597       }
598       ierr = ISDestroy(&is_subset_B);CHKERRQ(ierr);
599       ierr = MatCreateSchurComplement(AE_II,AE_II,AE_IE,AE_EI,AE_EE,&S_Ej);CHKERRQ(ierr);
600       ierr = MatDestroy(&AE_EE);CHKERRQ(ierr);
601       ierr = MatDestroy(&AE_IE);CHKERRQ(ierr);
602       ierr = MatDestroy(&AE_EI);CHKERRQ(ierr);
603       if (AE_II == A_II) { /* we can reuse the same ksp */
604         KSP ksp;
605         ierr = MatSchurComplementGetKSP(sub_schurs->S,&ksp);CHKERRQ(ierr);
606         ierr = MatSchurComplementSetKSP(S_Ej,ksp);CHKERRQ(ierr);
607       } else { /* build new ksp object which inherits ksp and pc types from the original one */
608         KSP       origksp,schurksp;
609         PC        origpc,schurpc;
610         KSPType   ksp_type;
611         PetscInt  n_internal;
612         PetscBool ispcnone;
613 
614         ierr = MatSchurComplementGetKSP(sub_schurs->S,&origksp);CHKERRQ(ierr);
615         ierr = MatSchurComplementGetKSP(S_Ej,&schurksp);CHKERRQ(ierr);
616         ierr = KSPGetType(origksp,&ksp_type);CHKERRQ(ierr);
617         ierr = KSPSetType(schurksp,ksp_type);CHKERRQ(ierr);
618         ierr = KSPGetPC(schurksp,&schurpc);CHKERRQ(ierr);
619         ierr = KSPGetPC(origksp,&origpc);CHKERRQ(ierr);
620         ierr = PetscObjectTypeCompare((PetscObject)origpc,PCNONE,&ispcnone);CHKERRQ(ierr);
621         if (!ispcnone) {
622           PCType pc_type;
623           ierr = PCGetType(origpc,&pc_type);CHKERRQ(ierr);
624           ierr = PCSetType(schurpc,pc_type);CHKERRQ(ierr);
625         } else {
626           ierr = PCSetType(schurpc,PCLU);CHKERRQ(ierr);
627         }
628         ierr = ISGetSize(is_I,&n_internal);CHKERRQ(ierr);
629         if (n_internal) { /* UMFPACK gives error with 0 sized problems */
630           MatSolverPackage solver=NULL;
631           ierr = PCFactorGetMatSolverPackage(origpc,(const MatSolverPackage*)&solver);CHKERRQ(ierr);
632           if (solver) {
633             ierr = PCFactorSetMatSolverPackage(schurpc,solver);CHKERRQ(ierr);
634           }
635         }
636         ierr = KSPSetUp(schurksp);CHKERRQ(ierr);
637       }
638       ierr = ISGetLocalSize(sub_schurs->is_subs[i],&subset_size);CHKERRQ(ierr);
639       ierr = MatCreateSeqDense(PETSC_COMM_SELF,subset_size,subset_size,work,&S_Ej_expl);CHKERRQ(ierr);
640       ierr = PCBDDCComputeExplicitSchur(S_Ej,sub_schurs->is_hermitian,MAT_REUSE_MATRIX,&S_Ej_expl);CHKERRQ(ierr);
641       ierr = PetscObjectTypeCompare((PetscObject)S_Ej_expl,MATSEQDENSE,&Sdense);CHKERRQ(ierr);
642       if (Sdense) {
643         for (j=0;j<subset_size;j++) {
644           dummy_idx[j]=local_size+j;
645         }
646         ierr = MatSetValues(sub_schurs->S_Ej_all,subset_size,dummy_idx,subset_size,dummy_idx,work,INSERT_VALUES);CHKERRQ(ierr);
647       } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not yet implemented for sparse matrices");
648       ierr = MatDestroy(&S_Ej);CHKERRQ(ierr);
649       ierr = MatDestroy(&S_Ej_expl);CHKERRQ(ierr);
650       local_size += subset_size;
651     }
652     ierr = PetscFree2(dummy_idx,work);CHKERRQ(ierr);
653     /* free */
654     ierr = ISDestroy(&is_I);CHKERRQ(ierr);
655     ierr = MatDestroy(&AE_II);CHKERRQ(ierr);
656     ierr = PetscFree(all_local_idx_N);CHKERRQ(ierr);
657   } else {
658     Mat         A;
659     IS          is_A_all;
660     PetscScalar *work,*S_data;
661     PetscInt    *idxs_schur,n_I,n_I_all,*dummy_idx,size_schur,size_active_schur,cum,cum2;
662     PetscBool   mumps_S;
663 
664     /* get working mat */
665     n_I = 0;
666     if (is_I_layer) {
667       ierr = ISGetLocalSize(is_I_layer,&n_I);CHKERRQ(ierr);
668     }
669     if (!sub_schurs->is_dir) {
670       ierr = ISCreateGeneral(PETSC_COMM_SELF,local_size+n_I,all_local_idx_N,PETSC_COPY_VALUES,&is_A_all);CHKERRQ(ierr);
671       size_schur = local_size;
672     } else {
673       IS list[2];
674 
675       ierr = ISCreateGeneral(PETSC_COMM_SELF,local_size+n_I,all_local_idx_N,PETSC_COPY_VALUES,&list[0]);CHKERRQ(ierr);
676       list[1] = sub_schurs->is_dir;
677       ierr = ISConcatenate(PETSC_COMM_SELF,2,list,&is_A_all);CHKERRQ(ierr);
678       ierr = ISDestroy(&list[0]);CHKERRQ(ierr);
679       ierr = ISGetLocalSize(sub_schurs->is_dir,&size_schur);CHKERRQ(ierr);
680       size_schur += local_size;
681     }
682     ierr = PetscFree(all_local_idx_N);CHKERRQ(ierr);
683     size_active_schur = local_size; /* size active schurs does not count any dirichlet dof on the interface */
684     ierr = MatGetSubMatrix(sub_schurs->A,is_A_all,is_A_all,MAT_INITIAL_MATRIX,&A);CHKERRQ(ierr);
685     ierr = MatSetOptionsPrefix(A,"sub_schurs_");CHKERRQ(ierr);
686     ierr = MatSetOption(A,MAT_SYMMETRIC,sub_schurs->is_hermitian);CHKERRQ(ierr);
687     ierr = MatSetOption(A,MAT_HERMITIAN,sub_schurs->is_hermitian);CHKERRQ(ierr);
688     ierr = MatSetOption(A,MAT_SPD,sub_schurs->is_posdef);CHKERRQ(ierr);
689 
690     if (n_I) {
691       if (sub_schurs->is_hermitian && sub_schurs->is_posdef) {
692         ierr = MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_CHOLESKY,&F);CHKERRQ(ierr);
693       } else {
694         ierr = MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_LU,&F);CHKERRQ(ierr);
695       }
696       /* subsets ordered last */
697       ierr = PetscMalloc1(size_schur,&idxs_schur);CHKERRQ(ierr);
698       for (i=0;i<size_schur;i++) {
699         idxs_schur[i] = n_I+i+1;
700       }
701 #if defined(PETSC_HAVE_MUMPS)
702       ierr = MatMumpsSetSchurIndices(F,size_schur,idxs_schur);CHKERRQ(ierr);
703 #endif
704       ierr = PetscFree(idxs_schur);CHKERRQ(ierr);
705 
706       /* factorization step */
707       if (sub_schurs->is_hermitian && sub_schurs->is_posdef) {
708         ierr = MatCholeskyFactorSymbolic(F,A,NULL,NULL);CHKERRQ(ierr);
709 #if defined(PETSC_HAVE_MUMPS) /* be sure that icntl 19 is not set by command line */
710         ierr = MatMumpsSetIcntl(F,19,2);CHKERRQ(ierr);
711 #endif
712         ierr = MatCholeskyFactorNumeric(F,A,NULL);CHKERRQ(ierr);
713       } else {
714         ierr = MatLUFactorSymbolic(F,A,NULL,NULL,NULL);CHKERRQ(ierr);
715 #if defined(PETSC_HAVE_MUMPS) /* be sure that icntl 19 is not set by command line */
716         ierr = MatMumpsSetIcntl(F,19,3);CHKERRQ(ierr);
717 #endif
718         ierr = MatLUFactorNumeric(F,A,NULL);CHKERRQ(ierr);
719       }
720 
721       /* get explicit Schur Complement computed during numeric factorization */
722 #if defined(PETSC_HAVE_MUMPS)
723       ierr = MatMumpsGetSchurComplement(F,&S_all);CHKERRQ(ierr);
724 #endif
725 
726       /* we can reuse the solvers if we are not using the economic version */
727       ierr = ISGetLocalSize(sub_schurs->is_I,&n_I_all);CHKERRQ(ierr);
728       reuse_solvers = (PetscBool)(reuse_solvers && (n_I == n_I_all));
729       mumps_S = PETSC_TRUE;
730     } else { /* we can't use MUMPS when size_schur == size_of_the_problem */
731       ierr = MatConvert(A,MATSEQDENSE,MAT_INITIAL_MATRIX,&S_all);CHKERRQ(ierr);
732       reuse_solvers = PETSC_FALSE;
733       mumps_S = PETSC_FALSE;
734     }
735 
736     if (reuse_solvers) {
737       Mat              A_II;
738       Vec              vec1_B;
739       PCBDDCReuseMumps msolv_ctx;
740 
741       if (sub_schurs->reuse_mumps) {
742         ierr = PCBDDCReuseMumpsReset(sub_schurs->reuse_mumps);CHKERRQ(ierr);
743       } else {
744         ierr = PetscNew(&sub_schurs->reuse_mumps);CHKERRQ(ierr);
745       }
746       msolv_ctx = sub_schurs->reuse_mumps;
747       ierr = MatSchurComplementGetSubMatrices(sub_schurs->S,&A_II,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
748       ierr = MatGetSize(A_II,&msolv_ctx->n_I,NULL);CHKERRQ(ierr);
749       ierr = PetscObjectReference((PetscObject)F);CHKERRQ(ierr);
750       msolv_ctx->F = F;
751       ierr = MatCreateVecs(F,&msolv_ctx->sol,&msolv_ctx->rhs);CHKERRQ(ierr);
752 
753       /* interior solver */
754       ierr = PCCreate(PETSC_COMM_SELF,&msolv_ctx->interior_solver);CHKERRQ(ierr);
755       ierr = PCSetOperators(msolv_ctx->interior_solver,A_II,A_II);CHKERRQ(ierr);
756       ierr = PCSetType(msolv_ctx->interior_solver,PCSHELL);CHKERRQ(ierr);
757       ierr = PCShellSetContext(msolv_ctx->interior_solver,msolv_ctx);CHKERRQ(ierr);
758       ierr = PCShellSetApply(msolv_ctx->interior_solver,PCBDDCMumpsInteriorSolve);CHKERRQ(ierr);
759       ierr = PCShellSetApplyTranspose(msolv_ctx->interior_solver,PCBDDCMumpsInteriorSolveTranspose);CHKERRQ(ierr);
760 
761       /* correction solver */
762       ierr = PCCreate(PETSC_COMM_SELF,&msolv_ctx->correction_solver);CHKERRQ(ierr);
763       ierr = PCSetOperators(msolv_ctx->correction_solver,A,A);CHKERRQ(ierr);
764       ierr = PCSetType(msolv_ctx->correction_solver,PCSHELL);CHKERRQ(ierr);
765       ierr = PCShellSetContext(msolv_ctx->correction_solver,msolv_ctx);CHKERRQ(ierr);
766       ierr = PCShellSetApply(msolv_ctx->correction_solver,PCBDDCMumpsCorrectionSolve);CHKERRQ(ierr);
767       ierr = PCShellSetApplyTranspose(msolv_ctx->correction_solver,PCBDDCMumpsCorrectionSolveTranspose);CHKERRQ(ierr);
768 
769       /* scatter and vecs for Schur complement solver */
770       ierr = MatCreateVecs(S_all,&msolv_ctx->sol_B,&msolv_ctx->rhs_B);CHKERRQ(ierr);
771       ierr = MatCreateVecs(sub_schurs->S,&vec1_B,NULL);CHKERRQ(ierr);
772       ierr = ISGlobalToLocalMappingApplyIS(sub_schurs->BtoNmap,IS_GTOLM_DROP,is_A_all,&msolv_ctx->is_B);CHKERRQ(ierr);
773       ierr = VecScatterCreate(vec1_B,msolv_ctx->is_B,msolv_ctx->sol_B,NULL,&msolv_ctx->correction_scatter_B);CHKERRQ(ierr);
774       ierr = VecDestroy(&vec1_B);CHKERRQ(ierr);
775       ierr = PetscObjectReference((PetscObject)is_A_all);CHKERRQ(ierr);
776       msolv_ctx->is_R = is_A_all;
777     }
778     ierr = MatDestroy(&A);CHKERRQ(ierr);
779     ierr = ISDestroy(&is_A_all);CHKERRQ(ierr);
780 
781     /* Work arrays */
782     ierr = PetscMalloc2(max_subset_size,&dummy_idx,max_subset_size*max_subset_size,&work);CHKERRQ(ierr);
783 
784     /* matrices for adaptive selection */
785     if (compute_Stilda) {
786       if (!sub_schurs->sum_S_Ej_tilda_all) {
787         ierr = MatCreate(PETSC_COMM_SELF,&sub_schurs->sum_S_Ej_tilda_all);CHKERRQ(ierr);
788         ierr = MatSetSizes(sub_schurs->sum_S_Ej_tilda_all,PETSC_DECIDE,PETSC_DECIDE,size_active_schur,size_active_schur);CHKERRQ(ierr);
789         ierr = MatSetType(sub_schurs->sum_S_Ej_tilda_all,MATAIJ);CHKERRQ(ierr);
790         ierr = MatSeqAIJSetPreallocation(sub_schurs->sum_S_Ej_tilda_all,0,nnz);CHKERRQ(ierr);
791       }
792       if (!sub_schurs->sum_S_Ej_inv_all) {
793         ierr = MatCreate(PETSC_COMM_SELF,&sub_schurs->sum_S_Ej_inv_all);CHKERRQ(ierr);
794         ierr = MatSetSizes(sub_schurs->sum_S_Ej_inv_all,PETSC_DECIDE,PETSC_DECIDE,size_active_schur,size_active_schur);CHKERRQ(ierr);
795         ierr = MatSetType(sub_schurs->sum_S_Ej_inv_all,MATAIJ);CHKERRQ(ierr);
796         ierr = MatSeqAIJSetPreallocation(sub_schurs->sum_S_Ej_inv_all,0,nnz);CHKERRQ(ierr);
797       }
798     }
799 
800     /* S_Ej_all */
801     cum = cum2 = 0;
802     ierr = MatDenseGetArray(S_all,&S_data);CHKERRQ(ierr);
803     for (i=0;i<sub_schurs->n_subs;i++) {
804       PetscInt j;
805 
806       /* get S_E */
807       ierr = ISGetLocalSize(sub_schurs->is_subs[i],&subset_size);CHKERRQ(ierr);
808       if (mumps_S && sub_schurs->is_hermitian) { /* When using MUMPS data I need to expand to upper triangular (column oriented) */
809         PetscInt k;
810         for (k=0;k<subset_size;k++) {
811           for (j=k;j<subset_size;j++) {
812             work[k*subset_size+j] = S_data[cum2+k*size_schur+j];
813             work[j*subset_size+k] = S_data[cum2+k*size_schur+j];
814           }
815         }
816       } else { /* copy to workspace */
817         PetscInt k;
818         for (k=0;k<subset_size;k++) {
819           for (j=0;j<subset_size;j++) {
820             work[k*subset_size+j] = S_data[cum2+k*size_schur+j];
821           }
822         }
823       }
824       /* insert S_E values */
825       for (j=0;j<subset_size;j++) dummy_idx[j] = cum+j;
826       ierr = MatSetValues(sub_schurs->S_Ej_all,subset_size,dummy_idx,subset_size,dummy_idx,work,INSERT_VALUES);CHKERRQ(ierr);
827 
828       /* if adaptivity is requested, invert S_E block */
829       if (compute_Stilda) {
830         ierr = PetscBLASIntCast(subset_size,&B_N);CHKERRQ(ierr);
831         ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
832         if (sub_schurs->is_hermitian && sub_schurs->is_posdef) { /* TODO add sytrf/i for symmetric non hermitian */
833           PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&B_N,work,&B_N,&B_ierr));
834           if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in POTRF Lapack routine %d",(int)B_ierr);
835           PetscStackCallBLAS("LAPACKpotri",LAPACKpotri_("L",&B_N,work,&B_N,&B_ierr));
836           if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in POTRI Lapack routine %d",(int)B_ierr);
837         } else {
838           PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&B_N,&B_N,work,&B_N,pivots,&B_ierr));
839           if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRF Lapack routine %d",(int)B_ierr);
840           PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&B_N,work,&B_N,pivots,Bwork,&B_lwork,&B_ierr));
841           if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRI Lapack routine %d",(int)B_ierr);
842         }
843         ierr = PetscFPTrapPop();CHKERRQ(ierr);
844         ierr = MatSetValues(sub_schurs->sum_S_Ej_inv_all,subset_size,dummy_idx,subset_size,dummy_idx,work,INSERT_VALUES);CHKERRQ(ierr);
845       }
846       cum += subset_size;
847       cum2 += subset_size*(size_schur + 1);
848     }
849     ierr = MatDenseRestoreArray(S_all,&S_data);CHKERRQ(ierr);
850 
851 #if defined(PETSC_HAVE_MUMPS)
852     if (mumps_S) {
853       ierr = MatMumpsRestoreSchurComplement(F,&S_all);CHKERRQ(ierr);
854     }
855 #endif
856 
857     if (compute_Stilda && size_active_schur) {
858       if (sub_schurs->n_subs == 1 && size_schur == size_active_schur) { /* we already computed the inverse */
859         PetscInt j;
860         for (j=0;j<size_schur;j++) dummy_idx[j] = j;
861         ierr = MatSetValues(sub_schurs->sum_S_Ej_tilda_all,size_schur,dummy_idx,size_schur,dummy_idx,work,INSERT_VALUES);CHKERRQ(ierr);
862       } else {
863         if (mumps_S) { /* use MatMumps calls to invert S */
864 #if defined(PETSC_HAVE_MUMPS)
865           ierr = MatMumpsInvertSchurComplement(F);CHKERRQ(ierr);
866           ierr = MatMumpsGetSchurComplement(F,&S_all);CHKERRQ(ierr);
867 #endif
868         } else { /* we need to invert explicitly since we are not using MUMPS for S */
869           ierr = MatDenseGetArray(S_all,&S_data);CHKERRQ(ierr);
870           ierr = PetscBLASIntCast(size_schur,&B_N);CHKERRQ(ierr);
871           ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
872           if (sub_schurs->is_hermitian && sub_schurs->is_posdef) { /* TODO add sytrf/i for symmetric non hermitian */
873             PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&B_N,S_data,&B_N,&B_ierr));
874             if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in POTRF Lapack routine %d",(int)B_ierr);
875             PetscStackCallBLAS("LAPACKpotri",LAPACKpotri_("L",&B_N,S_data,&B_N,&B_ierr));
876             if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in POTRI Lapack routine %d",(int)B_ierr);
877           } else {
878             PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&B_N,&B_N,S_data,&B_N,pivots,&B_ierr));
879             if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRF Lapack routine %d",(int)B_ierr);
880             PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&B_N,S_data,&B_N,pivots,Bwork,&B_lwork,&B_ierr));
881             if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRI Lapack routine %d",(int)B_ierr);
882           }
883           ierr = PetscFPTrapPop();CHKERRQ(ierr);
884           ierr = MatDenseRestoreArray(S_all,&S_data);CHKERRQ(ierr);
885         }
886         /* S_Ej_tilda_all */
887         cum = cum2 = 0;
888         ierr = MatDenseGetArray(S_all,&S_data);CHKERRQ(ierr);
889         for (i=0;i<sub_schurs->n_subs;i++) {
890           PetscInt j;
891 
892           ierr = ISGetLocalSize(sub_schurs->is_subs[i],&subset_size);CHKERRQ(ierr);
893           /* get (St^-1)_E */
894           if (sub_schurs->is_hermitian) { /* Here I don't need to expand to upper triangular (column oriented) */
895             PetscInt k;
896             for (k=0;k<subset_size;k++) {
897               for (j=k;j<subset_size;j++) {
898                 work[k*subset_size+j] = S_data[cum2+k*size_schur+j];
899               }
900             }
901           } else { /* copy to workspace */
902             PetscInt k;
903             for (k=0;k<subset_size;k++) {
904               for (j=0;j<subset_size;j++) {
905                 work[k*subset_size+j] = S_data[cum2+k*size_schur+j];
906               }
907             }
908           }
909           for (j=0;j<subset_size;j++) dummy_idx[j] = cum+j;
910           ierr = MatSetValues(sub_schurs->sum_S_Ej_tilda_all,subset_size,dummy_idx,subset_size,dummy_idx,work,INSERT_VALUES);CHKERRQ(ierr);
911           cum += subset_size;
912           cum2 += subset_size*(size_schur + 1);
913         }
914         ierr = MatDenseRestoreArray(S_all,&S_data);CHKERRQ(ierr);
915 #if defined(PETSC_HAVE_MUMPS)
916         if (mumps_S) {
917           ierr = MatMumpsRestoreSchurComplement(F,&S_all);CHKERRQ(ierr);
918         }
919 #endif
920       }
921     }
922     ierr = PetscFree2(dummy_idx,work);CHKERRQ(ierr);
923   }
924   ierr = PetscFree(nnz);CHKERRQ(ierr);
925   ierr = MatDestroy(&F);CHKERRQ(ierr);
926   ierr = ISDestroy(&is_I_layer);CHKERRQ(ierr);
927   ierr = MatDestroy(&S_all);CHKERRQ(ierr);
928   ierr = MatDestroy(&A_BB);CHKERRQ(ierr);
929   ierr = MatDestroy(&A_IB);CHKERRQ(ierr);
930   ierr = MatDestroy(&A_BI);CHKERRQ(ierr);
931   ierr = MatAssemblyBegin(sub_schurs->S_Ej_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
932   ierr = MatAssemblyEnd(sub_schurs->S_Ej_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
933   if (compute_Stilda) {
934     ierr = MatAssemblyBegin(sub_schurs->sum_S_Ej_tilda_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
935     ierr = MatAssemblyEnd(sub_schurs->sum_S_Ej_tilda_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
936     ierr = MatAssemblyBegin(sub_schurs->sum_S_Ej_inv_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
937     ierr = MatAssemblyEnd(sub_schurs->sum_S_Ej_inv_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
938   }
939 
940   /* Global matrix of all assembled Schur on subsets */
941   ierr = MatISSetLocalMat(work_mat,sub_schurs->S_Ej_all);CHKERRQ(ierr);
942   ierr = MatISSetMPIXAIJPreallocation_Private(work_mat,global_schur_subsets,PETSC_TRUE);CHKERRQ(ierr);
943   ierr = MatISGetMPIXAIJ(work_mat,MAT_REUSE_MATRIX,&global_schur_subsets);CHKERRQ(ierr);
944 
945   /* Get local part of (\sum_j S_Ej) */
946   if (!sub_schurs->sum_S_Ej_all) {
947     ierr = MatGetSubMatrices(global_schur_subsets,1,&all_subsets_n,&all_subsets_n,MAT_INITIAL_MATRIX,&submats);CHKERRQ(ierr);
948     sub_schurs->sum_S_Ej_all = submats[0];
949   } else {
950     ierr = PetscMalloc1(1,&submats);CHKERRQ(ierr);
951     submats[0] = sub_schurs->sum_S_Ej_all;
952     ierr = MatGetSubMatrices(global_schur_subsets,1,&all_subsets_n,&all_subsets_n,MAT_REUSE_MATRIX,&submats);CHKERRQ(ierr);
953   }
954 
955   /* Compute explicitly (\sum_j S_Ej)^-1 (faster scaling during PCApply, needs extra work when doing setup) */
956   if (faster_deluxe) {
957     Mat         tmpmat;
958     PetscScalar *array;
959     PetscInt    cum;
960 
961     ierr = MatSeqAIJGetArray(sub_schurs->sum_S_Ej_all,&array);CHKERRQ(ierr);
962     cum = 0;
963     for (i=0;i<sub_schurs->n_subs;i++) {
964       ierr = ISGetLocalSize(sub_schurs->is_subs[i],&subset_size);CHKERRQ(ierr);
965       ierr = PetscBLASIntCast(subset_size,&B_N);CHKERRQ(ierr);
966       ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
967       if (sub_schurs->is_hermitian && sub_schurs->is_posdef) {
968         PetscInt j,k;
969 
970         PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&B_N,array+cum,&B_N,&B_ierr));
971         if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in POTRF Lapack routine %d",(int)B_ierr);
972         PetscStackCallBLAS("LAPACKpotri",LAPACKpotri_("L",&B_N,array+cum,&B_N,&B_ierr));
973         if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in POTRI Lapack routine %d",(int)B_ierr);
974         for (j=0;j<B_N;j++) {
975           for (k=j+1;k<B_N;k++) {
976             array[k*B_N+j+cum] = array[j*B_N+k+cum];
977           }
978         }
979       } else {
980         PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&B_N,&B_N,array+cum,&B_N,pivots,&B_ierr));
981         if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRF Lapack routine %d",(int)B_ierr);
982         PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&B_N,array+cum,&B_N,pivots,Bwork,&B_lwork,&B_ierr));
983         if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRI Lapack routine %d",(int)B_ierr);
984       }
985       ierr = PetscFPTrapPop();CHKERRQ(ierr);
986       cum += subset_size*subset_size;
987     }
988     ierr = MatSeqAIJRestoreArray(sub_schurs->sum_S_Ej_all,&array);CHKERRQ(ierr);
989     ierr = MatMatMult(sub_schurs->S_Ej_all,sub_schurs->sum_S_Ej_all,MAT_INITIAL_MATRIX,1.0,&tmpmat);CHKERRQ(ierr);
990     ierr = MatDestroy(&sub_schurs->S_Ej_all);CHKERRQ(ierr);
991     ierr = MatDestroy(&sub_schurs->sum_S_Ej_all);CHKERRQ(ierr);
992     sub_schurs->S_Ej_all = tmpmat;
993   }
994 
995   /* Get local part of (\sum_j S^-1_Ej) (\sum_j St^-1_Ej) */
996   if (compute_Stilda) {
997     ierr = MatISSetLocalMat(work_mat,sub_schurs->sum_S_Ej_tilda_all);CHKERRQ(ierr);
998     ierr = MatISGetMPIXAIJ(work_mat,MAT_REUSE_MATRIX,&global_schur_subsets);CHKERRQ(ierr);
999     submats[0] = sub_schurs->sum_S_Ej_tilda_all;
1000     ierr = MatGetSubMatrices(global_schur_subsets,1,&all_subsets_n,&all_subsets_n,MAT_REUSE_MATRIX,&submats);CHKERRQ(ierr);
1001     ierr = MatISSetLocalMat(work_mat,sub_schurs->sum_S_Ej_inv_all);CHKERRQ(ierr);
1002     ierr = MatISGetMPIXAIJ(work_mat,MAT_REUSE_MATRIX,&global_schur_subsets);CHKERRQ(ierr);
1003     submats[0] = sub_schurs->sum_S_Ej_inv_all;
1004     ierr = MatGetSubMatrices(global_schur_subsets,1,&all_subsets_n,&all_subsets_n,MAT_REUSE_MATRIX,&submats);CHKERRQ(ierr);
1005   }
1006 
1007   /* free workspace */
1008   ierr = PetscFree(submats);CHKERRQ(ierr);
1009   ierr = PetscFree2(Bwork,pivots);CHKERRQ(ierr);
1010   ierr = MatDestroy(&global_schur_subsets);CHKERRQ(ierr);
1011   ierr = MatDestroy(&work_mat);CHKERRQ(ierr);
1012   ierr = ISDestroy(&all_subsets_n);CHKERRQ(ierr);
1013   ierr = PetscCommDestroy(&comm_n);CHKERRQ(ierr);
1014   PetscFunctionReturn(0);
1015 }
1016 
1017 #undef __FUNCT__
1018 #define __FUNCT__ "PCBDDCSubSchursInit"
1019 PetscErrorCode PCBDDCSubSchursInit(PCBDDCSubSchurs sub_schurs, IS is_I, IS is_B, PCBDDCGraph graph, ISLocalToGlobalMapping BtoNmap)
1020 {
1021   IS              *faces,*edges,*all_cc,vertices;
1022   PetscInt        i,n_faces,n_edges,n_all_cc;
1023   PetscBool       is_sorted;
1024   PetscErrorCode  ierr;
1025 
1026   PetscFunctionBegin;
1027   ierr = ISSorted(is_I,&is_sorted);CHKERRQ(ierr);
1028   if (!is_sorted) SETERRQ(PetscObjectComm((PetscObject)is_I),PETSC_ERR_PLIB,"IS for I dofs should be shorted");
1029   ierr = ISSorted(is_B,&is_sorted);CHKERRQ(ierr);
1030   if (!is_sorted) SETERRQ(PetscObjectComm((PetscObject)is_B),PETSC_ERR_PLIB,"IS for B dofs should be shorted");
1031 
1032   /* reset any previous data */
1033   ierr = PCBDDCSubSchursReset(sub_schurs);CHKERRQ(ierr);
1034 
1035   /* get index sets for faces and edges (already sorted by global ordering) */
1036   ierr = PCBDDCGraphGetCandidatesIS(graph,&n_faces,&faces,&n_edges,&edges,&vertices);CHKERRQ(ierr);
1037   n_all_cc = n_faces+n_edges;
1038   ierr = PetscBTCreate(n_all_cc,&sub_schurs->is_edge);CHKERRQ(ierr);
1039   ierr = PetscMalloc1(n_all_cc,&all_cc);CHKERRQ(ierr);
1040   for (i=0;i<n_faces;i++) {
1041     all_cc[i] = faces[i];
1042   }
1043   for (i=0;i<n_edges;i++) {
1044     all_cc[n_faces+i] = edges[i];
1045     ierr = PetscBTSet(sub_schurs->is_edge,n_faces+i);CHKERRQ(ierr);
1046   }
1047   ierr = PetscFree(faces);CHKERRQ(ierr);
1048   ierr = PetscFree(edges);CHKERRQ(ierr);
1049   sub_schurs->is_dir = NULL;
1050   ierr = PCBDDCGraphGetDirichletDofsB(graph,&sub_schurs->is_dir);CHKERRQ(ierr);
1051 
1052   /* Determine if MUMPS can be used */
1053   sub_schurs->use_mumps = PETSC_FALSE;
1054 #if defined(PETSC_HAVE_MUMPS)
1055   sub_schurs->use_mumps = PETSC_TRUE;
1056 #endif
1057 
1058   ierr = PetscObjectReference((PetscObject)is_I);CHKERRQ(ierr);
1059   sub_schurs->is_I = is_I;
1060   ierr = PetscObjectReference((PetscObject)is_B);CHKERRQ(ierr);
1061   sub_schurs->is_B = is_B;
1062   ierr = PetscObjectReference((PetscObject)graph->l2gmap);CHKERRQ(ierr);
1063   sub_schurs->l2gmap = graph->l2gmap;
1064   ierr = PetscObjectReference((PetscObject)BtoNmap);CHKERRQ(ierr);
1065   sub_schurs->BtoNmap = BtoNmap;
1066   sub_schurs->n_subs = n_all_cc;
1067   sub_schurs->is_subs = all_cc;
1068   if (!sub_schurs->use_mumps) { /* sort by local ordering mumps is not present */
1069     for (i=0;i<sub_schurs->n_subs;i++) {
1070       ierr = ISSort(sub_schurs->is_subs[i]);CHKERRQ(ierr);
1071     }
1072   }
1073   sub_schurs->is_vertices = vertices;
1074   sub_schurs->S_Ej_all = NULL;
1075   sub_schurs->sum_S_Ej_all = NULL;
1076   sub_schurs->sum_S_Ej_inv_all = NULL;
1077   sub_schurs->sum_S_Ej_tilda_all = NULL;
1078   sub_schurs->is_Ej_all = NULL;
1079   PetscFunctionReturn(0);
1080 }
1081 
1082 #undef __FUNCT__
1083 #define __FUNCT__ "PCBDDCSubSchursCreate"
1084 PetscErrorCode PCBDDCSubSchursCreate(PCBDDCSubSchurs *sub_schurs)
1085 {
1086   PCBDDCSubSchurs schurs_ctx;
1087   PetscErrorCode  ierr;
1088 
1089   PetscFunctionBegin;
1090   ierr = PetscNew(&schurs_ctx);CHKERRQ(ierr);
1091   schurs_ctx->n_subs = 0;
1092   *sub_schurs = schurs_ctx;
1093   PetscFunctionReturn(0);
1094 }
1095 
1096 #undef __FUNCT__
1097 #define __FUNCT__ "PCBDDCSubSchursDestroy"
1098 PetscErrorCode PCBDDCSubSchursDestroy(PCBDDCSubSchurs *sub_schurs)
1099 {
1100   PetscErrorCode ierr;
1101 
1102   PetscFunctionBegin;
1103   ierr = PCBDDCSubSchursReset(*sub_schurs);CHKERRQ(ierr);
1104   ierr = PetscFree(*sub_schurs);CHKERRQ(ierr);
1105   PetscFunctionReturn(0);
1106 }
1107 
1108 #undef __FUNCT__
1109 #define __FUNCT__ "PCBDDCSubSchursReset"
1110 PetscErrorCode PCBDDCSubSchursReset(PCBDDCSubSchurs sub_schurs)
1111 {
1112   PetscInt       i;
1113   PetscErrorCode ierr;
1114 
1115   PetscFunctionBegin;
1116   ierr = MatDestroy(&sub_schurs->A);CHKERRQ(ierr);
1117   ierr = MatDestroy(&sub_schurs->S);CHKERRQ(ierr);
1118   ierr = ISDestroy(&sub_schurs->is_I);CHKERRQ(ierr);
1119   ierr = ISDestroy(&sub_schurs->is_B);CHKERRQ(ierr);
1120   ierr = ISLocalToGlobalMappingDestroy(&sub_schurs->l2gmap);CHKERRQ(ierr);
1121   ierr = ISLocalToGlobalMappingDestroy(&sub_schurs->BtoNmap);CHKERRQ(ierr);
1122   ierr = MatDestroy(&sub_schurs->S_Ej_all);CHKERRQ(ierr);
1123   ierr = MatDestroy(&sub_schurs->sum_S_Ej_all);CHKERRQ(ierr);
1124   ierr = MatDestroy(&sub_schurs->sum_S_Ej_inv_all);CHKERRQ(ierr);
1125   ierr = MatDestroy(&sub_schurs->sum_S_Ej_tilda_all);CHKERRQ(ierr);
1126   ierr = ISDestroy(&sub_schurs->is_Ej_all);CHKERRQ(ierr);
1127   ierr = ISDestroy(&sub_schurs->is_vertices);CHKERRQ(ierr);
1128   ierr = ISDestroy(&sub_schurs->is_dir);CHKERRQ(ierr);
1129   ierr = PetscBTDestroy(&sub_schurs->is_edge);CHKERRQ(ierr);
1130   for (i=0;i<sub_schurs->n_subs;i++) {
1131     ierr = ISDestroy(&sub_schurs->is_subs[i]);CHKERRQ(ierr);
1132   }
1133   if (sub_schurs->n_subs) {
1134     ierr = PetscFree(sub_schurs->is_subs);CHKERRQ(ierr);
1135   }
1136   if (sub_schurs->reuse_mumps) {
1137     ierr = PCBDDCReuseMumpsReset(sub_schurs->reuse_mumps);CHKERRQ(ierr);
1138   }
1139   ierr = PetscFree(sub_schurs->reuse_mumps);CHKERRQ(ierr);
1140   sub_schurs->n_subs = 0;
1141   PetscFunctionReturn(0);
1142 }
1143 
1144 #undef __FUNCT__
1145 #define __FUNCT__ "PCBDDCAdjGetNextLayer_Private"
1146 PETSC_STATIC_INLINE PetscErrorCode PCBDDCAdjGetNextLayer_Private(PetscInt* queue_tip,PetscInt n_prev,PetscBT touched,PetscInt* xadj,PetscInt* adjncy,PetscInt* n_added)
1147 {
1148   PetscInt       i,j,n;
1149   PetscErrorCode ierr;
1150 
1151   PetscFunctionBegin;
1152   n = 0;
1153   for (i=-n_prev;i<0;i++) {
1154     PetscInt start_dof = queue_tip[i];
1155     for (j=xadj[start_dof];j<xadj[start_dof+1];j++) {
1156       PetscInt dof = adjncy[j];
1157       if (!PetscBTLookup(touched,dof)) {
1158         ierr = PetscBTSet(touched,dof);CHKERRQ(ierr);
1159         queue_tip[n] = dof;
1160         n++;
1161       }
1162     }
1163   }
1164   *n_added = n;
1165   PetscFunctionReturn(0);
1166 }
1167