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