xref: /petsc/src/ksp/pc/impls/bddc/bddcschurs.c (revision b1b3d7a2dc467bc01d0ce7eb62bfacc6b4178a73) 
1 #include <../src/ksp/pc/impls/bddc/bddc.h>
2 #include <../src/ksp/pc/impls/bddc/bddcprivate.h>
3 
4 static PetscErrorCode PCBDDCAdjGetNextLayer_Private(PetscInt*,PetscInt,PetscBT,PetscInt*,PetscInt*,PetscInt*);
5 
6 #undef __FUNCT__
7 #define __FUNCT__ "PCBDDCSubSchursSetUpNew"
8 PetscErrorCode PCBDDCSubSchursSetUpNew(PCBDDCSubSchurs sub_schurs, PetscInt xadj[], PetscInt adjncy[], PetscInt nlayers)
9 {
10   Mat                    A_II,A_IB,A_BI,A_BB;
11   Mat                    AE_II,*AE_IE,*AE_EI,*AE_EE;
12   IS                     is_I,*is_subset_B;
13   ISLocalToGlobalMapping BtoNmap;
14   PetscInt               i;
15   PetscBool              implicit_schurs;
16   PetscErrorCode         ierr;
17 
18   PetscFunctionBegin;
19   ierr = PetscObjectTypeCompare((PetscObject)sub_schurs->S,MATSCHURCOMPLEMENT,&implicit_schurs);CHKERRQ(ierr);
20 
21   /* allocate space for schur complements */
22   ierr = PetscMalloc5(sub_schurs->n_subs,&sub_schurs->is_AEj_I,
23                       sub_schurs->n_subs,&sub_schurs->is_AEj_B,
24                       sub_schurs->n_subs,&sub_schurs->S_Ej,
25                       sub_schurs->n_subs,&sub_schurs->work1,
26                       sub_schurs->n_subs,&sub_schurs->work2);CHKERRQ(ierr);
27 
28   /* get Schur complement matrices */
29   if (implicit_schurs) {
30     ierr = MatSchurComplementGetSubMatrices(sub_schurs->S,&A_II,NULL,&A_IB,&A_BI,&A_BB);CHKERRQ(ierr);
31     ierr = PetscMalloc4(sub_schurs->n_subs,&is_subset_B,
32                         sub_schurs->n_subs,&AE_IE,
33                         sub_schurs->n_subs,&AE_EI,
34                         sub_schurs->n_subs,&AE_EE);CHKERRQ(ierr);
35   }
36 
37   /* determine interior problems */
38   if (nlayers >= 0 && xadj != NULL && adjncy != NULL) { /* Interior problems can be different from the original one */
39     PetscBT                touched;
40     const PetscInt*        idx_B;
41     PetscInt               n_I,n_B,n_local_dofs,n_prev_added,j,layer,*local_numbering;
42 
43     /* get sizes */
44     ierr = ISGetLocalSize(sub_schurs->is_I,&n_I);CHKERRQ(ierr);
45     ierr = ISGetLocalSize(sub_schurs->is_B,&n_B);CHKERRQ(ierr);
46 
47     ierr = PetscMalloc1(n_I+n_B,&local_numbering);CHKERRQ(ierr);
48     ierr = PetscBTCreate(n_I+n_B,&touched);CHKERRQ(ierr);
49     ierr = PetscBTMemzero(n_I+n_B,touched);CHKERRQ(ierr);
50 
51     /* all boundary dofs must be skipped when adding layers */
52     ierr = ISGetIndices(sub_schurs->is_B,&idx_B);CHKERRQ(ierr);
53     for (j=0;j<n_B;j++) {
54       ierr = PetscBTSet(touched,idx_B[j]);CHKERRQ(ierr);
55     }
56     ierr = PetscMemcpy(local_numbering,idx_B,n_B*sizeof(PetscInt));CHKERRQ(ierr);
57     ierr = ISRestoreIndices(sub_schurs->is_B,&idx_B);CHKERRQ(ierr);
58 
59     /* add prescribed number of layers of dofs */
60     n_local_dofs = n_B;
61     n_prev_added = n_B;
62     for (layer=0;layer<nlayers;layer++) {
63       PetscInt n_added;
64       if (n_local_dofs == n_I+n_B) break;
65       if (n_local_dofs > n_I+n_B) {
66         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);
67       }
68       ierr = PCBDDCAdjGetNextLayer_Private(local_numbering+n_local_dofs,n_prev_added,touched,xadj,adjncy,&n_added);CHKERRQ(ierr);
69       n_prev_added = n_added;
70       n_local_dofs += n_added;
71       if (!n_added) break;
72     }
73     ierr = PetscBTDestroy(&touched);CHKERRQ(ierr);
74 
75     /* IS for I dofs in original numbering */
76     ierr = ISCreateGeneral(PetscObjectComm((PetscObject)sub_schurs->is_I),n_local_dofs-n_B,local_numbering+n_B,PETSC_COPY_VALUES,&sub_schurs->is_AEj_I[0]);CHKERRQ(ierr);
77     ierr = PetscFree(local_numbering);CHKERRQ(ierr);
78     ierr = ISSort(sub_schurs->is_AEj_I[0]);CHKERRQ(ierr);
79     /* IS for I dofs in boundary numbering */
80     if (implicit_schurs) {
81       ISLocalToGlobalMapping ItoNmap;
82       ierr = ISLocalToGlobalMappingCreateIS(sub_schurs->is_I,&ItoNmap);CHKERRQ(ierr);
83       ierr = ISGlobalToLocalMappingApplyIS(ItoNmap,IS_GTOLM_DROP,sub_schurs->is_AEj_I[0],&is_I);CHKERRQ(ierr);
84       ierr = ISLocalToGlobalMappingDestroy(&ItoNmap);CHKERRQ(ierr);
85 
86       /* II block */
87       ierr = MatGetSubMatrix(A_II,is_I,is_I,MAT_INITIAL_MATRIX,&AE_II);CHKERRQ(ierr);
88     }
89   } else {
90     PetscInt n_I;
91 
92     /* IS for I dofs in original numbering */
93     ierr = PetscObjectReference((PetscObject)sub_schurs->is_I);CHKERRQ(ierr);
94     sub_schurs->is_AEj_I[0] = sub_schurs->is_I;
95 
96     /* IS for I dofs in I numbering (strided 1) */
97     if (implicit_schurs) {
98       ierr = ISGetSize(sub_schurs->is_I,&n_I);CHKERRQ(ierr);
99       ierr = ISCreateStride(PetscObjectComm((PetscObject)sub_schurs->is_I),n_I,0,1,&is_I);CHKERRQ(ierr);
100 
101       /* II block is the same */
102       ierr = PetscObjectReference((PetscObject)A_II);CHKERRQ(ierr);
103       AE_II = A_II;
104     }
105   }
106 
107   /* TODO: just for compatibility with the previous version, needs to be fixed */
108   for (i=1;i<sub_schurs->n_subs;i++) {
109     ierr = PetscObjectReference((PetscObject)sub_schurs->is_AEj_I[0]);CHKERRQ(ierr);
110     sub_schurs->is_AEj_I[i] = sub_schurs->is_AEj_I[0];
111   }
112 
113   if (implicit_schurs) {
114     /* subsets in original and boundary numbering */
115     ierr = ISLocalToGlobalMappingCreateIS(sub_schurs->is_B,&BtoNmap);CHKERRQ(ierr);
116     for (i=0;i<sub_schurs->n_subs;i++) {
117       ierr = ISDuplicate(sub_schurs->is_subs[i],&sub_schurs->is_AEj_B[i]);CHKERRQ(ierr);
118       ierr = ISSort(sub_schurs->is_AEj_B[i]);CHKERRQ(ierr);
119       ierr = ISGlobalToLocalMappingApplyIS(BtoNmap,IS_GTOLM_DROP,sub_schurs->is_AEj_B[i],&is_subset_B[i]);CHKERRQ(ierr);
120     }
121     ierr = ISLocalToGlobalMappingDestroy(&BtoNmap);CHKERRQ(ierr);
122 
123     /* EE block */
124     for (i=0;i<sub_schurs->n_subs;i++) {
125       ierr = MatGetSubMatrix(A_BB,is_subset_B[i],is_subset_B[i],MAT_INITIAL_MATRIX,&AE_EE[i]);CHKERRQ(ierr);
126     }
127     /* IE block */
128     for (i=0;i<sub_schurs->n_subs;i++) {
129       ierr = MatGetSubMatrix(A_IB,is_I,is_subset_B[i],MAT_INITIAL_MATRIX,&AE_IE[i]);CHKERRQ(ierr);
130     }
131     /* EI block */
132     for (i=0;i<sub_schurs->n_subs;i++) {
133       ierr = MatGetSubMatrix(A_BI,is_subset_B[i],is_I,MAT_INITIAL_MATRIX,&AE_EI[i]);CHKERRQ(ierr);
134     }
135 
136     /* setup Schur complements on subset */
137     for (i=0;i<sub_schurs->n_subs;i++) {
138       ierr = MatCreateSchurComplement(AE_II,AE_II,AE_IE[i],AE_EI[i],AE_EE[i],&sub_schurs->S_Ej[i]);CHKERRQ(ierr);
139       ierr = MatCreateVecs(sub_schurs->S_Ej[i],&sub_schurs->work1[i],&sub_schurs->work2[i]);CHKERRQ(ierr);
140       if (AE_II == A_II) { /* we can reuse the same ksp */
141         KSP ksp;
142         ierr = MatSchurComplementGetKSP(sub_schurs->S,&ksp);CHKERRQ(ierr);
143         ierr = MatSchurComplementSetKSP(sub_schurs->S_Ej[i],ksp);CHKERRQ(ierr);
144       } else { /* build new ksp object which inherits ksp and pc types from the original one */
145         KSP      origksp,schurksp;
146         PC       origpc,schurpc;
147         KSPType  ksp_type;
148         PCType   pc_type;
149         PetscInt n_internal;
150 
151         ierr = MatSchurComplementGetKSP(sub_schurs->S,&origksp);CHKERRQ(ierr);
152         ierr = MatSchurComplementGetKSP(sub_schurs->S_Ej[i],&schurksp);CHKERRQ(ierr);
153         ierr = KSPGetType(origksp,&ksp_type);CHKERRQ(ierr);
154         ierr = KSPSetType(schurksp,ksp_type);CHKERRQ(ierr);
155         ierr = KSPGetPC(schurksp,&schurpc);CHKERRQ(ierr);
156         ierr = KSPGetPC(origksp,&origpc);CHKERRQ(ierr);
157         ierr = PCGetType(origpc,&pc_type);CHKERRQ(ierr);
158         ierr = PCSetType(schurpc,pc_type);CHKERRQ(ierr);
159         ierr = ISGetSize(is_I,&n_internal);CHKERRQ(ierr);
160         if (n_internal) { /* UMFPACK gives error with 0 sized problems */
161           MatSolverPackage solver=NULL;
162           ierr = PCFactorGetMatSolverPackage(origpc,(const MatSolverPackage*)&solver);CHKERRQ(ierr);
163           if (solver) {
164             ierr = PCFactorSetMatSolverPackage(schurpc,solver);CHKERRQ(ierr);
165           }
166         }
167         ierr = KSPSetUp(schurksp);CHKERRQ(ierr);
168       }
169     }
170     /* free */
171     ierr = ISDestroy(&is_I);CHKERRQ(ierr);
172     ierr = MatDestroy(&AE_II);CHKERRQ(ierr);
173     for (i=0;i<sub_schurs->n_subs;i++) {
174       ierr = MatDestroy(&AE_EE[i]);CHKERRQ(ierr);
175       ierr = MatDestroy(&AE_IE[i]);CHKERRQ(ierr);
176       ierr = MatDestroy(&AE_EI[i]);CHKERRQ(ierr);
177       ierr = ISDestroy(&is_subset_B[i]);CHKERRQ(ierr);
178     }
179     ierr = PetscFree4(is_subset_B,AE_IE,AE_EI,AE_EE);CHKERRQ(ierr);
180   }
181   PetscFunctionReturn(0);
182 }
183 
184 #undef __FUNCT__
185 #define __FUNCT__ "PCBDDCSubSchursInit"
186 PetscErrorCode PCBDDCSubSchursInit(PCBDDCSubSchurs sub_schurs, Mat S, IS is_I, IS is_B, PCBDDCGraph graph, PetscInt seqthreshold)
187 {
188   IS                  *faces,*edges,*all_cc;
189   PetscInt            *index_sequential,*index_parallel;
190   PetscInt            *auxlocal_sequential,*auxlocal_parallel;
191   PetscInt            *auxglobal_sequential,*auxglobal_parallel;
192   PetscInt            *auxmapping;//,*idxs;
193   PetscInt            i,max_subset_size;
194   PetscInt            n_sequential_problems,n_local_sequential_problems,n_parallel_problems,n_local_parallel_problems;
195   PetscInt            n_faces,n_edges,n_all_cc;
196   PetscBool           is_sorted;
197   PetscErrorCode  ierr;
198 
199   PetscFunctionBegin;
200   ierr = ISSorted(is_I,&is_sorted);CHKERRQ(ierr);
201   if (!is_sorted) {
202     SETERRQ(PetscObjectComm((PetscObject)is_I),PETSC_ERR_PLIB,"IS for I dofs should be shorted");
203   }
204   ierr = ISSorted(is_B,&is_sorted);CHKERRQ(ierr);
205   if (!is_sorted) {
206     SETERRQ(PetscObjectComm((PetscObject)is_B),PETSC_ERR_PLIB,"IS for B dofs should be shorted");
207   }
208 
209   /* reset any previous data */
210   ierr = PCBDDCSubSchursReset(sub_schurs);CHKERRQ(ierr);
211 
212   /* get index sets for faces and edges */
213   ierr = PCBDDCGraphGetCandidatesIS(graph,&n_faces,&faces,&n_edges,&edges,NULL);CHKERRQ(ierr);
214   n_all_cc = n_faces+n_edges;
215   ierr = PetscMalloc1(n_all_cc,&all_cc);CHKERRQ(ierr);
216   for (i=0;i<n_faces;i++) {
217     all_cc[i] = faces[i];
218   }
219   for (i=0;i<n_edges;i++) {
220     all_cc[n_faces+i] = edges[i];
221   }
222   ierr = PetscFree(faces);CHKERRQ(ierr);
223   ierr = PetscFree(edges);CHKERRQ(ierr);
224 
225   /* map interface's subsets */
226   max_subset_size = 0;
227   for (i=0;i<n_all_cc;i++) {
228     PetscInt subset_size;
229     ierr = ISGetLocalSize(all_cc[i],&subset_size);CHKERRQ(ierr);
230     max_subset_size = PetscMax(max_subset_size,subset_size);
231   }
232   ierr = PetscMalloc1(max_subset_size,&auxmapping);CHKERRQ(ierr);
233   ierr = PetscMalloc2(graph->ncc,&auxlocal_sequential,
234                       graph->ncc,&auxlocal_parallel);CHKERRQ(ierr);
235   ierr = PetscMalloc2(graph->ncc,&index_sequential,
236                       graph->ncc,&index_parallel);CHKERRQ(ierr);
237 
238   /* if threshold is negative, uses all sequential problems */
239   if (seqthreshold < 0) seqthreshold = max_subset_size;
240 
241   /* determine which problem has to be solved in parallel or sequentially */
242   n_local_sequential_problems = 0;
243   n_local_parallel_problems = 0;
244   for (i=0;i<n_all_cc;i++) {
245     PetscInt       subset_size,j,min_loc = 0;
246     const PetscInt *idxs;
247 
248     ierr = ISGetLocalSize(all_cc[i],&subset_size);CHKERRQ(ierr);
249     ierr = ISGetIndices(all_cc[i],&idxs);CHKERRQ(ierr);
250     ierr = ISLocalToGlobalMappingApply(graph->l2gmap,subset_size,idxs,auxmapping);CHKERRQ(ierr);
251     for (j=1;j<subset_size;j++) {
252       if (auxmapping[j]<auxmapping[min_loc]) {
253         min_loc = j;
254       }
255     }
256     if (subset_size > seqthreshold) {
257       index_parallel[n_local_parallel_problems] = i;
258       auxlocal_parallel[n_local_parallel_problems] = idxs[min_loc];
259       n_local_parallel_problems++;
260     } else {
261       index_sequential[n_local_sequential_problems] = i;
262       auxlocal_sequential[n_local_sequential_problems] = idxs[min_loc];
263       n_local_sequential_problems++;
264     }
265     ierr = ISRestoreIndices(all_cc[i],&idxs);CHKERRQ(ierr);
266   }
267 
268   /* Number parallel problems */
269   auxglobal_parallel = 0;
270   ierr = PCBDDCSubsetNumbering(PetscObjectComm((PetscObject)graph->l2gmap),graph->l2gmap,n_local_parallel_problems,auxlocal_parallel,PETSC_NULL,&n_parallel_problems,&auxglobal_parallel);CHKERRQ(ierr);
271 
272   /* Number sequential problems */
273   auxglobal_sequential = 0;
274   ierr = PCBDDCSubsetNumbering(PetscObjectComm((PetscObject)graph->l2gmap),graph->l2gmap,n_local_sequential_problems,auxlocal_sequential,PETSC_NULL,&n_sequential_problems,&auxglobal_sequential);CHKERRQ(ierr);
275 
276   /* update info in sub_schurs */
277   ierr = PetscObjectReference((PetscObject)S);CHKERRQ(ierr);
278   sub_schurs->S = S;
279   ierr = PetscObjectReference((PetscObject)is_I);CHKERRQ(ierr);
280   sub_schurs->is_I = is_I;
281   ierr = PetscObjectReference((PetscObject)is_B);CHKERRQ(ierr);
282   sub_schurs->is_B = is_B;
283   sub_schurs->n_subs_seq = n_local_sequential_problems;
284   sub_schurs->n_subs_par = n_local_parallel_problems;
285   sub_schurs->n_subs_seq_g = n_sequential_problems;
286   sub_schurs->n_subs_par_g = n_parallel_problems;
287   sub_schurs->n_subs = sub_schurs->n_subs_seq + sub_schurs->n_subs_par;
288   sub_schurs->is_subs = all_cc;
289   sub_schurs->index_sequential = index_sequential;
290   sub_schurs->index_parallel = index_parallel;
291   sub_schurs->auxglobal_sequential = auxglobal_sequential;
292   sub_schurs->auxglobal_parallel = auxglobal_parallel;
293 
294   /* free workspace */
295   ierr = PetscFree(auxmapping);CHKERRQ(ierr);
296   ierr = PetscFree2(auxlocal_sequential,auxlocal_parallel);CHKERRQ(ierr);
297 
298   PetscFunctionReturn(0);
299 }
300 
301 #undef __FUNCT__
302 #define __FUNCT__ "PCBDDCSubSchursCreate"
303 PetscErrorCode PCBDDCSubSchursCreate(PCBDDCSubSchurs *sub_schurs)
304 {
305   PCBDDCSubSchurs schurs_ctx;
306   PetscErrorCode  ierr;
307 
308   PetscFunctionBegin;
309   ierr = PetscNew(&schurs_ctx);CHKERRQ(ierr);
310   schurs_ctx->n_subs = 0;
311   *sub_schurs = schurs_ctx;
312   PetscFunctionReturn(0);
313 }
314 
315 #undef __FUNCT__
316 #define __FUNCT__ "PCBDDCSubSchursDestroy"
317 PetscErrorCode PCBDDCSubSchursDestroy(PCBDDCSubSchurs *sub_schurs)
318 {
319   PetscErrorCode ierr;
320 
321   PetscFunctionBegin;
322   ierr = PCBDDCSubSchursReset(*sub_schurs);CHKERRQ(ierr);
323   ierr = PetscFree(*sub_schurs);CHKERRQ(ierr);
324   PetscFunctionReturn(0);
325 }
326 
327 #undef __FUNCT__
328 #define __FUNCT__ "PCBDDCSubSchursReset"
329 PetscErrorCode PCBDDCSubSchursReset(PCBDDCSubSchurs sub_schurs)
330 {
331   PetscInt       i;
332   PetscErrorCode ierr;
333 
334   PetscFunctionBegin;
335   ierr = MatDestroy(&sub_schurs->S);CHKERRQ(ierr);
336   ierr = ISDestroy(&sub_schurs->is_I);CHKERRQ(ierr);
337   ierr = ISDestroy(&sub_schurs->is_B);CHKERRQ(ierr);
338   for (i=0;i<sub_schurs->n_subs;i++) {
339     ierr = ISDestroy(&sub_schurs->is_subs[i]);CHKERRQ(ierr);
340     ierr = ISDestroy(&sub_schurs->is_AEj_I[i]);CHKERRQ(ierr);
341     ierr = ISDestroy(&sub_schurs->is_AEj_B[i]);CHKERRQ(ierr);
342     ierr = MatDestroy(&sub_schurs->S_Ej[i]);CHKERRQ(ierr);
343     ierr = VecDestroy(&sub_schurs->work1[i]);CHKERRQ(ierr);
344     ierr = VecDestroy(&sub_schurs->work2[i]);CHKERRQ(ierr);
345   }
346   if (sub_schurs->n_subs) {
347     ierr = PetscFree(sub_schurs->is_subs);CHKERRQ(ierr);
348     ierr = PetscFree5(sub_schurs->is_AEj_I,sub_schurs->is_AEj_B,sub_schurs->S_Ej,sub_schurs->work1,sub_schurs->work2);CHKERRQ(ierr);
349     ierr = PetscFree2(sub_schurs->index_sequential,sub_schurs->index_parallel);CHKERRQ(ierr);
350     ierr = PetscFree(sub_schurs->auxglobal_sequential);CHKERRQ(ierr);
351     ierr = PetscFree(sub_schurs->auxglobal_parallel);CHKERRQ(ierr);
352   }
353   sub_schurs->n_subs = 0;
354   PetscFunctionReturn(0);
355 }
356 
357 #undef __FUNCT__
358 #define __FUNCT__ "PCBDDCAdjGetNextLayer_Private"
359 PETSC_STATIC_INLINE PetscErrorCode PCBDDCAdjGetNextLayer_Private(PetscInt* queue_tip,PetscInt n_prev,PetscBT touched,PetscInt* xadj,PetscInt* adjncy,PetscInt* n_added)
360 {
361   PetscInt       i,j,n;
362   PetscErrorCode ierr;
363 
364   PetscFunctionBegin;
365   n = 0;
366   for (i=-n_prev;i<0;i++) {
367     PetscInt start_dof = queue_tip[i];
368     for (j=xadj[start_dof];j<xadj[start_dof+1];j++) {
369       PetscInt dof = adjncy[j];
370       if (!PetscBTLookup(touched,dof)) {
371         ierr = PetscBTSet(touched,dof);CHKERRQ(ierr);
372         queue_tip[n] = dof;
373         n++;
374       }
375     }
376   }
377   *n_added = n;
378   PetscFunctionReturn(0);
379 }
380 
381 #undef __FUNCT__
382 #define __FUNCT__ "PCBDDCSubSchursSetUp"
383 PetscErrorCode PCBDDCSubSchursSetUp(PCBDDCSubSchurs sub_schurs, Mat S, IS is_A_I, IS is_A_B, PetscInt ncc, IS is_cc[], PetscInt xadj[], PetscInt adjncy[], PetscInt nlayers)
384 {
385   Mat                    A_II,A_IB,A_BI,A_BB;
386   Mat                    AE_II,*AE_IE,*AE_EI,*AE_EE;
387   IS                     is_I,*is_subset_B;
388   ISLocalToGlobalMapping BtoNmap;
389   PetscInt               i;
390   PetscBool              is_sorted;
391   PetscErrorCode         ierr;
392 
393   PetscFunctionBegin;
394   ierr = ISSorted(is_A_I,&is_sorted);CHKERRQ(ierr);
395   if (!is_sorted) {
396     SETERRQ(PetscObjectComm((PetscObject)is_A_I),PETSC_ERR_PLIB,"IS for I dofs should be shorted");
397   }
398   ierr = ISSorted(is_A_B,&is_sorted);CHKERRQ(ierr);
399   if (!is_sorted) {
400     SETERRQ(PetscObjectComm((PetscObject)is_A_B),PETSC_ERR_PLIB,"IS for B dofs should be shorted");
401   }
402 
403   /* get Schur complement matrices */
404   ierr = MatSchurComplementGetSubMatrices(S,&A_II,NULL,&A_IB,&A_BI,&A_BB);CHKERRQ(ierr);
405 
406   /* allocate space for schur complements */
407   ierr = PetscMalloc5(ncc,&sub_schurs->is_AEj_I,ncc,&sub_schurs->is_AEj_B,ncc,&sub_schurs->S_Ej,ncc,&sub_schurs->work1,ncc,&sub_schurs->work2);CHKERRQ(ierr);
408   ierr = PetscMalloc4(ncc,&is_subset_B,ncc,&AE_IE,ncc,&AE_EI,ncc,&AE_EE);CHKERRQ(ierr);
409   sub_schurs->n_subs = ncc;
410 
411   /* maps */
412   if (sub_schurs->n_subs && nlayers >= 0 && xadj != NULL && adjncy != NULL) { /* Interior problems can be different from the original one */
413     ISLocalToGlobalMapping ItoNmap;
414     PetscBT                touched;
415     const PetscInt*        idx_B;
416     PetscInt               n_I,n_B,n_local_dofs,n_prev_added,j,layer,*local_numbering;
417 
418     /* get sizes */
419     ierr = ISGetLocalSize(is_A_I,&n_I);CHKERRQ(ierr);
420     ierr = ISGetLocalSize(is_A_B,&n_B);CHKERRQ(ierr);
421 
422     ierr = ISLocalToGlobalMappingCreateIS(is_A_I,&ItoNmap);CHKERRQ(ierr);
423     ierr = PetscMalloc1(n_I+n_B,&local_numbering);CHKERRQ(ierr);
424     ierr = PetscBTCreate(n_I+n_B,&touched);CHKERRQ(ierr);
425     ierr = PetscBTMemzero(n_I+n_B,touched);CHKERRQ(ierr);
426 
427     /* all boundary dofs must be skipped when adding layers */
428     ierr = ISGetIndices(is_A_B,&idx_B);CHKERRQ(ierr);
429     for (j=0;j<n_B;j++) {
430       ierr = PetscBTSet(touched,idx_B[j]);CHKERRQ(ierr);
431     }
432     ierr = PetscMemcpy(local_numbering,idx_B,n_B*sizeof(PetscInt));CHKERRQ(ierr);
433     ierr = ISRestoreIndices(is_A_B,&idx_B);CHKERRQ(ierr);
434 
435     /* add next layers of dofs */
436     n_local_dofs = n_B;
437     n_prev_added = n_B;
438     for (layer=0;layer<nlayers;layer++) {
439       PetscInt n_added;
440       if (n_local_dofs == n_I+n_B) break;
441       if (n_local_dofs > n_I+n_B) {
442         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);
443       }
444       ierr = PCBDDCAdjGetNextLayer_Private(local_numbering+n_local_dofs,n_prev_added,touched,xadj,adjncy,&n_added);CHKERRQ(ierr);
445       n_prev_added = n_added;
446       n_local_dofs += n_added;
447       if (!n_added) break;
448     }
449     ierr = PetscBTDestroy(&touched);CHKERRQ(ierr);
450 
451     /* IS for I dofs in original numbering and in I numbering */
452     ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ItoNmap),n_local_dofs-n_B,local_numbering+n_B,PETSC_COPY_VALUES,&sub_schurs->is_AEj_I[0]);CHKERRQ(ierr);
453     ierr = PetscFree(local_numbering);CHKERRQ(ierr);
454     ierr = ISSort(sub_schurs->is_AEj_I[0]);CHKERRQ(ierr);
455     ierr = ISGlobalToLocalMappingApplyIS(ItoNmap,IS_GTOLM_DROP,sub_schurs->is_AEj_I[0],&is_I);CHKERRQ(ierr);
456     ierr = ISLocalToGlobalMappingDestroy(&ItoNmap);CHKERRQ(ierr);
457 
458     /* II block */
459     ierr = MatGetSubMatrix(A_II,is_I,is_I,MAT_INITIAL_MATRIX,&AE_II);CHKERRQ(ierr);
460   } else {
461     PetscInt n_I;
462 
463     /* IS for I dofs in original numbering */
464     ierr = PetscObjectReference((PetscObject)is_A_I);CHKERRQ(ierr);
465     sub_schurs->is_AEj_I[0] = is_A_I;
466 
467     /* IS for I dofs in I numbering (strided 1) */
468     ierr = ISGetSize(is_A_I,&n_I);CHKERRQ(ierr);
469     ierr = ISCreateStride(PetscObjectComm((PetscObject)is_A_I),n_I,0,1,&is_I);CHKERRQ(ierr);
470 
471     /* II block is the same */
472     ierr = PetscObjectReference((PetscObject)A_II);CHKERRQ(ierr);
473     AE_II = A_II;
474   }
475 
476   /* TODO: just for compatibility with the previous version, needs to be fixed */
477   for (i=1;i<sub_schurs->n_subs;i++) {
478     ierr = PetscObjectReference((PetscObject)sub_schurs->is_AEj_I[0]);CHKERRQ(ierr);
479     sub_schurs->is_AEj_I[i] = sub_schurs->is_AEj_I[0];
480   }
481 
482   /* subsets in original and boundary numbering */
483   ierr = ISLocalToGlobalMappingCreateIS(is_A_B,&BtoNmap);CHKERRQ(ierr);
484   for (i=0;i<sub_schurs->n_subs;i++) {
485     ierr = ISDuplicate(is_cc[i],&sub_schurs->is_AEj_B[i]);CHKERRQ(ierr);
486     ierr = ISSort(sub_schurs->is_AEj_B[i]);CHKERRQ(ierr);
487     ierr = ISGlobalToLocalMappingApplyIS(BtoNmap,IS_GTOLM_DROP,sub_schurs->is_AEj_B[i],&is_subset_B[i]);CHKERRQ(ierr);
488   }
489   ierr = ISLocalToGlobalMappingDestroy(&BtoNmap);CHKERRQ(ierr);
490 
491   /* EE block */
492   for (i=0;i<sub_schurs->n_subs;i++) {
493     ierr = MatGetSubMatrix(A_BB,is_subset_B[i],is_subset_B[i],MAT_INITIAL_MATRIX,&AE_EE[i]);CHKERRQ(ierr);
494   }
495   /* IE block */
496   for (i=0;i<sub_schurs->n_subs;i++) {
497     ierr = MatGetSubMatrix(A_IB,is_I,is_subset_B[i],MAT_INITIAL_MATRIX,&AE_IE[i]);CHKERRQ(ierr);
498   }
499   /* EI block */
500   for (i=0;i<sub_schurs->n_subs;i++) {
501     ierr = MatGetSubMatrix(A_BI,is_subset_B[i],is_I,MAT_INITIAL_MATRIX,&AE_EI[i]);CHKERRQ(ierr);
502   }
503 
504   /* setup Schur complements on subset */
505   for (i=0;i<sub_schurs->n_subs;i++) {
506     ierr = MatCreateSchurComplement(AE_II,AE_II,AE_IE[i],AE_EI[i],AE_EE[i],&sub_schurs->S_Ej[i]);CHKERRQ(ierr);
507     ierr = MatCreateVecs(sub_schurs->S_Ej[i],&sub_schurs->work1[i],&sub_schurs->work2[i]);CHKERRQ(ierr);
508     if (AE_II == A_II) { /* we can reuse the same ksp */
509       KSP ksp;
510       ierr = MatSchurComplementGetKSP(S,&ksp);CHKERRQ(ierr);
511       ierr = MatSchurComplementSetKSP(sub_schurs->S_Ej[i],ksp);CHKERRQ(ierr);
512     } else { /* build new ksp object which inherits ksp and pc types from the original one */
513       KSP      origksp,schurksp;
514       PC       origpc,schurpc;
515       KSPType  ksp_type;
516       PCType   pc_type;
517       PetscInt n_internal;
518 
519       ierr = MatSchurComplementGetKSP(S,&origksp);CHKERRQ(ierr);
520       ierr = MatSchurComplementGetKSP(sub_schurs->S_Ej[i],&schurksp);CHKERRQ(ierr);
521       ierr = KSPGetType(origksp,&ksp_type);CHKERRQ(ierr);
522       ierr = KSPSetType(schurksp,ksp_type);CHKERRQ(ierr);
523       ierr = KSPGetPC(schurksp,&schurpc);CHKERRQ(ierr);
524       ierr = KSPGetPC(origksp,&origpc);CHKERRQ(ierr);
525       ierr = PCGetType(origpc,&pc_type);CHKERRQ(ierr);
526       ierr = PCSetType(schurpc,pc_type);CHKERRQ(ierr);
527       ierr = ISGetSize(is_I,&n_internal);CHKERRQ(ierr);
528       if (n_internal) { /* UMFPACK gives error with 0 sized problems */
529         MatSolverPackage solver=NULL;
530         ierr = PCFactorGetMatSolverPackage(origpc,(const MatSolverPackage*)&solver);CHKERRQ(ierr);
531         if (solver) {
532           ierr = PCFactorSetMatSolverPackage(schurpc,solver);CHKERRQ(ierr);
533         }
534       }
535       ierr = KSPSetUp(schurksp);CHKERRQ(ierr);
536     }
537   }
538   /* free */
539   ierr = ISDestroy(&is_I);CHKERRQ(ierr);
540   ierr = MatDestroy(&AE_II);CHKERRQ(ierr);
541   for (i=0;i<sub_schurs->n_subs;i++) {
542     ierr = MatDestroy(&AE_EE[i]);CHKERRQ(ierr);
543     ierr = MatDestroy(&AE_IE[i]);CHKERRQ(ierr);
544     ierr = MatDestroy(&AE_EI[i]);CHKERRQ(ierr);
545     ierr = ISDestroy(&is_subset_B[i]);CHKERRQ(ierr);
546   }
547   ierr = PetscFree4(is_subset_B,AE_IE,AE_EI,AE_EE);CHKERRQ(ierr);
548   PetscFunctionReturn(0);
549 }
550