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