xref: /petsc/src/ksp/pc/impls/bddc/bddcschurs.c (revision 9bb4a8ca6688cafd8c6250e0011bc4b45f2808a0) 
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,"Computation of Stildas 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     /* Loop on local problems end compute Schur complements explicitly */
369     local_size = 0;
370     for (i=0;i<sub_schurs->n_subs_seq;i++) {
371       Mat       S_Ej_expl;
372       PetscInt  j,lpi = sub_schurs->index_sequential[i];
373       PetscBool Sdense;
374 
375       ierr = PCBDDCComputeExplicitSchur(sub_schurs->S_Ej[lpi],&S_Ej_expl);CHKERRQ(ierr);
376       ierr = ISGetLocalSize(sub_schurs->is_subs[lpi],&subset_size);CHKERRQ(ierr);
377       ierr = PetscObjectTypeCompare((PetscObject)S_Ej_expl,MATSEQDENSE,&Sdense);CHKERRQ(ierr);
378       if (Sdense) {
379         for (j=0;j<subset_size;j++) {
380           dummy_idx[j]=local_size+j;
381         }
382         ierr = MatDenseGetArray(S_Ej_expl,&fill_vals);CHKERRQ(ierr);
383         ierr = MatSetValues(sub_schurs->S_Ej_all,subset_size,dummy_idx,subset_size,dummy_idx,fill_vals,INSERT_VALUES);CHKERRQ(ierr);
384         ierr = MatDenseRestoreArray(S_Ej_expl,&fill_vals);CHKERRQ(ierr);
385       } else {
386         SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not yet implemented for sparse matrices");
387       }
388       local_size += subset_size;
389       ierr = MatDestroy(&sub_schurs->S_Ej[lpi]);CHKERRQ(ierr);
390       sub_schurs->S_Ej[lpi] = S_Ej_expl;
391     }
392     /* Stildas are not computed without mumps */
393     for (i=0;i<sub_schurs->n_subs;i++) {
394       ierr = MatDestroy(&sub_schurs->S_Ej_tilda[i]);CHKERRQ(ierr);
395     }
396     ierr = PetscFree(dummy_idx);CHKERRQ(ierr);
397   } else {
398     PetscInt  *dummy_idx,n_all;
399     PetscBool is_symmetric;
400 
401     if (compute_Stilda) {
402       ierr = MatCreate(PETSC_COMM_SELF,&S_Ej_tilda_all);CHKERRQ(ierr);
403       ierr = MatSetSizes(S_Ej_tilda_all,PETSC_DECIDE,PETSC_DECIDE,local_size,local_size);CHKERRQ(ierr);
404       ierr = MatSetType(S_Ej_tilda_all,MATAIJ);CHKERRQ(ierr);
405       ierr = MatSeqAIJSetPreallocation(S_Ej_tilda_all,0,nnz);CHKERRQ(ierr);
406     }
407     /* Work arrays */
408     ierr = PetscMalloc1(max_subset_size,&dummy_idx);CHKERRQ(ierr);
409 
410     /* Work arrays */
411     ierr = MatGetSize(S_all,&n_all,NULL);CHKERRQ(ierr);
412     local_size = 0;
413     ierr = MatIsSymmetric(sub_schurs->A,0.0,&is_symmetric);CHKERRQ(ierr);
414     for (i=0;i<sub_schurs->n_subs;i++) {
415       IS  is_E;
416       PetscScalar *vals;
417       PetscInt j;
418 
419       ierr = ISGetLocalSize(sub_schurs->is_subs[i],&subset_size);CHKERRQ(ierr);
420       for (j=0;j<subset_size;j++) {
421         dummy_idx[j]=local_size+j;
422       }
423       ierr = ISCreateStride(PetscObjectComm((PetscObject)sub_schurs->is_I),subset_size,local_size,1,&is_E);CHKERRQ(ierr);
424       if (sub_schurs->S_Ej[i]) {
425         ierr = MatGetSubMatrix(S_all,is_E,is_E,MAT_REUSE_MATRIX,&sub_schurs->S_Ej[i]);CHKERRQ(ierr);
426       } else {
427         ierr = MatGetSubMatrix(S_all,is_E,is_E,MAT_INITIAL_MATRIX,&sub_schurs->S_Ej[i]);CHKERRQ(ierr);
428       }
429       ierr = MatDenseGetArray(sub_schurs->S_Ej[i],&vals);CHKERRQ(ierr);
430       ierr = MatSetValues(sub_schurs->S_Ej_all,subset_size,dummy_idx,subset_size,dummy_idx,vals,INSERT_VALUES);CHKERRQ(ierr);
431       ierr = MatDenseRestoreArray(sub_schurs->S_Ej[i],&vals);CHKERRQ(ierr);
432 
433       if (compute_Stilda) {
434         Mat Stilda;
435         if (sub_schurs->n_subs == 1) {
436           ierr = PetscObjectReference((PetscObject)sub_schurs->S_Ej[i]);CHKERRQ(ierr);
437           Stilda = sub_schurs->S_Ej[i];
438         } else {
439           KSP ksp;
440           PC  pc;
441           Mat S_EF,S_FE,S_FF,Stilda_impl;
442           IS  is_F;
443           PetscScalar eps=1.e-8;
444           PetscBool chop=PETSC_FALSE;
445 
446           ierr = ISComplement(is_E,0,n_all,&is_F);CHKERRQ(ierr);
447           ierr = MatGetSubMatrix(S_all,is_E,is_F,MAT_INITIAL_MATRIX,&S_EF);CHKERRQ(ierr);
448           ierr = MatGetSubMatrix(S_all,is_F,is_F,MAT_INITIAL_MATRIX,&S_FF);CHKERRQ(ierr);
449           ierr = MatGetSubMatrix(S_all,is_F,is_E,MAT_INITIAL_MATRIX,&S_FE);CHKERRQ(ierr);
450           ierr = ISDestroy(&is_F);CHKERRQ(ierr);
451           if (chop) {
452             ierr = MatChop(S_FF,eps);CHKERRQ(ierr);
453             ierr = MatConvert(S_FF,MATAIJ,MAT_REUSE_MATRIX,&S_FF);CHKERRQ(ierr);
454           }
455           ierr = MatCreateSchurComplement(S_FF,S_FF,S_FE,S_EF,sub_schurs->S_Ej[i],&Stilda_impl);CHKERRQ(ierr);
456           ierr = MatDestroy(&S_FF);CHKERRQ(ierr);
457           ierr = MatDestroy(&S_FE);CHKERRQ(ierr);
458           ierr = MatDestroy(&S_EF);CHKERRQ(ierr);
459           ierr = MatSchurComplementGetKSP(Stilda_impl,&ksp);CHKERRQ(ierr);
460           ierr = KSPSetType(ksp,KSPPREONLY);CHKERRQ(ierr);
461           ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
462           if (is_symmetric) {
463             ierr = PCSetType(pc,PCCHOLESKY);CHKERRQ(ierr);
464           } else {
465             ierr = PCSetType(pc,PCLU);CHKERRQ(ierr);
466           }
467           if (!chop) {
468             ierr = PCFactorSetUseInPlace(pc,PETSC_TRUE);CHKERRQ(ierr);
469           } else {
470             ierr = PCFactorSetMatSolverPackage(pc,MATSOLVERMUMPS);CHKERRQ(ierr);
471           }
472           ierr = KSPSetUp(ksp);CHKERRQ(ierr);
473           ierr = PCBDDCComputeExplicitSchur(Stilda_impl,&Stilda);CHKERRQ(ierr);
474           ierr = MatDestroy(&Stilda_impl);CHKERRQ(ierr);
475         }
476 /*
477         PetscViewerPushFormat(PETSC_VIEWER_STDOUT_SELF,PETSC_VIEWER_ASCII_MATLAB);
478         ierr = MatView(Stilda,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
479         PetscViewerPopFormat(PETSC_VIEWER_STDOUT_SELF);
480 */
481         ierr = MatDestroy(&sub_schurs->S_Ej_tilda[i]);CHKERRQ(ierr);
482         sub_schurs->S_Ej_tilda[i] = Stilda;
483         ierr = MatDenseGetArray(sub_schurs->S_Ej_tilda[i],&vals);CHKERRQ(ierr);
484         ierr = MatSetValues(S_Ej_tilda_all,subset_size,dummy_idx,subset_size,dummy_idx,vals,INSERT_VALUES);CHKERRQ(ierr);
485         ierr = MatDenseRestoreArray(sub_schurs->S_Ej_tilda[i],&vals);CHKERRQ(ierr);
486       } else {
487         ierr = MatDestroy(&sub_schurs->S_Ej_tilda[i]);CHKERRQ(ierr);
488       }
489       ierr = ISDestroy(&is_E);CHKERRQ(ierr);
490       local_size += subset_size;
491     }
492     ierr = PetscFree(dummy_idx);CHKERRQ(ierr);
493   }
494   ierr = PetscFree(nnz);CHKERRQ(ierr);
495   ierr = MatDestroy(&S_all);CHKERRQ(ierr);
496   ierr = MatAssemblyBegin(sub_schurs->S_Ej_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
497   ierr = MatAssemblyEnd(sub_schurs->S_Ej_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
498   if (compute_Stilda) {
499     ierr = MatAssemblyBegin(S_Ej_tilda_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
500     ierr = MatAssemblyEnd(S_Ej_tilda_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
501   }
502 
503   /* Global matrix of all assembled Schur on subsets */
504   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);
505   ierr = ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)sub_schurs->l2gmap),1,local_size,all_local_idx_G,PETSC_COPY_VALUES,&l2gmap_subsets);CHKERRQ(ierr);
506   ierr = MatCreateIS(PetscObjectComm((PetscObject)sub_schurs->l2gmap),1,PETSC_DECIDE,PETSC_DECIDE,global_size,global_size,l2gmap_subsets,&work_mat);CHKERRQ(ierr);
507   ierr = ISLocalToGlobalMappingDestroy(&l2gmap_subsets);CHKERRQ(ierr);
508   ierr = MatISSetLocalMat(work_mat,sub_schurs->S_Ej_all);CHKERRQ(ierr);
509   ierr = MatISGetMPIXAIJ(work_mat,MAT_INITIAL_MATRIX,&global_schur_subsets);CHKERRQ(ierr);
510   /* Get local part of (\sum_j S_Ej) */
511   ierr = PetscFree(all_local_idx_N);CHKERRQ(ierr);
512   ierr = ISCreateGeneral(PetscObjectComm((PetscObject)sub_schurs->l2gmap),local_size,all_local_idx_G,PETSC_OWN_POINTER,&temp_is);CHKERRQ(ierr);
513   if (sub_schurs->sum_S_Ej_all) {
514     ierr = MatGetSubMatrixUnsorted(global_schur_subsets,temp_is,temp_is,MAT_REUSE_MATRIX,&sub_schurs->sum_S_Ej_all);CHKERRQ(ierr);
515   } else {
516     ierr = MatGetSubMatrixUnsorted(global_schur_subsets,temp_is,temp_is,MAT_INITIAL_MATRIX,&sub_schurs->sum_S_Ej_all);CHKERRQ(ierr);
517   }
518   if (compute_Stilda) {
519     ierr = MatISSetLocalMat(work_mat,S_Ej_tilda_all);CHKERRQ(ierr);
520     ierr = MatISGetMPIXAIJ(work_mat,MAT_REUSE_MATRIX,&global_schur_subsets);CHKERRQ(ierr);
521     if (sub_schurs->sum_S_Ej_tilda_all) {
522       ierr = MatGetSubMatrixUnsorted(global_schur_subsets,temp_is,temp_is,MAT_REUSE_MATRIX,&sub_schurs->sum_S_Ej_tilda_all);CHKERRQ(ierr);
523     } else {
524       ierr = MatGetSubMatrixUnsorted(global_schur_subsets,temp_is,temp_is,MAT_INITIAL_MATRIX,&sub_schurs->sum_S_Ej_tilda_all);CHKERRQ(ierr);
525     }
526   }
527 
528   ierr = MatDestroy(&global_schur_subsets);CHKERRQ(ierr);
529   ierr = MatDestroy(&S_Ej_tilda_all);CHKERRQ(ierr);
530   ierr = MatDestroy(&work_mat);CHKERRQ(ierr);
531   ierr = ISDestroy(&temp_is);CHKERRQ(ierr);
532   PetscFunctionReturn(0);
533 }
534 
535 #undef __FUNCT__
536 #define __FUNCT__ "PCBDDCSubSchursInit"
537 PetscErrorCode PCBDDCSubSchursInit(PCBDDCSubSchurs sub_schurs, Mat A, Mat S, IS is_I, IS is_B, PCBDDCGraph graph, ISLocalToGlobalMapping BtoNmap, PetscInt seqthreshold)
538 {
539   IS                  *faces,*edges,*all_cc,vertices;
540   PetscInt            *index_sequential,*index_parallel;
541   PetscInt            *auxlocal_sequential,*auxlocal_parallel;
542   PetscInt            *auxglobal_sequential,*auxglobal_parallel;
543   PetscInt            *auxmapping;
544   PetscInt            i,max_subset_size;
545   PetscInt            n_sequential_problems,n_local_sequential_problems,n_parallel_problems,n_local_parallel_problems;
546   PetscInt            n_faces,n_edges,n_all_cc;
547   PetscBool           is_sorted;
548   PetscErrorCode  ierr;
549 
550   PetscFunctionBegin;
551   ierr = ISSorted(is_I,&is_sorted);CHKERRQ(ierr);
552   if (!is_sorted) {
553     SETERRQ(PetscObjectComm((PetscObject)is_I),PETSC_ERR_PLIB,"IS for I dofs should be shorted");
554   }
555   ierr = ISSorted(is_B,&is_sorted);CHKERRQ(ierr);
556   if (!is_sorted) {
557     SETERRQ(PetscObjectComm((PetscObject)is_B),PETSC_ERR_PLIB,"IS for B dofs should be shorted");
558   }
559 
560   /* reset any previous data */
561   ierr = PCBDDCSubSchursReset(sub_schurs);CHKERRQ(ierr);
562 
563   /* get index sets for faces and edges */
564   ierr = PCBDDCGraphGetCandidatesIS(graph,&n_faces,&faces,&n_edges,&edges,&vertices);CHKERRQ(ierr);
565   n_all_cc = n_faces+n_edges;
566   ierr = PetscMalloc1(n_all_cc,&all_cc);CHKERRQ(ierr);
567   for (i=0;i<n_faces;i++) {
568     all_cc[i] = faces[i];
569   }
570   for (i=0;i<n_edges;i++) {
571     all_cc[n_faces+i] = edges[i];
572   }
573   ierr = PetscFree(faces);CHKERRQ(ierr);
574   ierr = PetscFree(edges);CHKERRQ(ierr);
575 
576   /* map interface's subsets */
577   max_subset_size = 0;
578   for (i=0;i<n_all_cc;i++) {
579     PetscInt subset_size;
580     ierr = ISGetLocalSize(all_cc[i],&subset_size);CHKERRQ(ierr);
581     max_subset_size = PetscMax(max_subset_size,subset_size);
582   }
583   ierr = PetscMalloc1(max_subset_size,&auxmapping);CHKERRQ(ierr);
584   ierr = PetscMalloc2(graph->ncc,&auxlocal_sequential,
585                       graph->ncc,&auxlocal_parallel);CHKERRQ(ierr);
586   ierr = PetscMalloc2(graph->ncc,&index_sequential,
587                       graph->ncc,&index_parallel);CHKERRQ(ierr);
588 
589   /* if threshold is negative uses all sequential problems (possibly using MUMPS) */
590   sub_schurs->use_mumps = PETSC_FALSE;
591   if (seqthreshold < 0) {
592     seqthreshold = max_subset_size;
593 #if defined(PETSC_HAVE_MUMPS)
594     sub_schurs->use_mumps = !!A;
595 #endif
596   }
597 
598 
599   /* determine which problem has to be solved in parallel or sequentially */
600   n_local_sequential_problems = 0;
601   n_local_parallel_problems = 0;
602   for (i=0;i<n_all_cc;i++) {
603     PetscInt       subset_size,j,min_loc = 0;
604     const PetscInt *idxs;
605 
606     ierr = ISGetLocalSize(all_cc[i],&subset_size);CHKERRQ(ierr);
607     ierr = ISGetIndices(all_cc[i],&idxs);CHKERRQ(ierr);
608     ierr = ISLocalToGlobalMappingApply(graph->l2gmap,subset_size,idxs,auxmapping);CHKERRQ(ierr);
609     for (j=1;j<subset_size;j++) {
610       if (auxmapping[j]<auxmapping[min_loc]) {
611         min_loc = j;
612       }
613     }
614     if (subset_size > seqthreshold) {
615       index_parallel[n_local_parallel_problems] = i;
616       auxlocal_parallel[n_local_parallel_problems] = idxs[min_loc];
617       n_local_parallel_problems++;
618     } else {
619       index_sequential[n_local_sequential_problems] = i;
620       auxlocal_sequential[n_local_sequential_problems] = idxs[min_loc];
621       n_local_sequential_problems++;
622     }
623     ierr = ISRestoreIndices(all_cc[i],&idxs);CHKERRQ(ierr);
624   }
625 
626   /* Number parallel problems */
627   auxglobal_parallel = 0;
628   ierr = PCBDDCSubsetNumbering(PetscObjectComm((PetscObject)graph->l2gmap),graph->l2gmap,n_local_parallel_problems,auxlocal_parallel,PETSC_NULL,&n_parallel_problems,&auxglobal_parallel);CHKERRQ(ierr);
629 
630   /* Number sequential problems */
631   auxglobal_sequential = 0;
632   ierr = PCBDDCSubsetNumbering(PetscObjectComm((PetscObject)graph->l2gmap),graph->l2gmap,n_local_sequential_problems,auxlocal_sequential,PETSC_NULL,&n_sequential_problems,&auxglobal_sequential);CHKERRQ(ierr);
633 
634   /* update info in sub_schurs */
635   if (sub_schurs->use_mumps && A) {
636     ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr);
637     sub_schurs->A = A;
638   }
639   ierr = PetscObjectReference((PetscObject)S);CHKERRQ(ierr);
640   sub_schurs->S = S;
641   ierr = PetscObjectReference((PetscObject)is_I);CHKERRQ(ierr);
642   sub_schurs->is_I = is_I;
643   ierr = PetscObjectReference((PetscObject)is_B);CHKERRQ(ierr);
644   sub_schurs->is_B = is_B;
645   ierr = PetscObjectReference((PetscObject)graph->l2gmap);CHKERRQ(ierr);
646   sub_schurs->l2gmap = graph->l2gmap;
647   ierr = PetscObjectReference((PetscObject)BtoNmap);CHKERRQ(ierr);
648   sub_schurs->BtoNmap = BtoNmap;
649   sub_schurs->n_subs_seq = n_local_sequential_problems;
650   sub_schurs->n_subs_par = n_local_parallel_problems;
651   sub_schurs->n_subs_seq_g = n_sequential_problems;
652   sub_schurs->n_subs_par_g = n_parallel_problems;
653   sub_schurs->n_subs = sub_schurs->n_subs_seq + sub_schurs->n_subs_par;
654   sub_schurs->is_subs = all_cc;
655   if (!sub_schurs->use_mumps) { /* for adaptive selection */
656     for (i=0;i<sub_schurs->n_subs;i++) {
657       ierr = ISSort(sub_schurs->is_subs[i]);CHKERRQ(ierr);
658     }
659   }
660   sub_schurs->is_Ej_com = vertices;
661   sub_schurs->index_sequential = index_sequential;
662   sub_schurs->index_parallel = index_parallel;
663   sub_schurs->auxglobal_sequential = auxglobal_sequential;
664   sub_schurs->auxglobal_parallel = auxglobal_parallel;
665 
666 
667   /* allocate space for schur complements */
668   ierr = PetscCalloc2(sub_schurs->n_subs,&sub_schurs->S_Ej,
669                       sub_schurs->n_subs,&sub_schurs->S_Ej_tilda);CHKERRQ(ierr);
670   sub_schurs->S_Ej_all = NULL;
671   sub_schurs->sum_S_Ej_all = NULL;
672   sub_schurs->sum_S_Ej_tilda_all = NULL;
673   sub_schurs->is_Ej_all = NULL;
674 
675   /* free workspace */
676   ierr = PetscFree(auxmapping);CHKERRQ(ierr);
677   ierr = PetscFree2(auxlocal_sequential,auxlocal_parallel);CHKERRQ(ierr);
678 
679   PetscFunctionReturn(0);
680 }
681 
682 #undef __FUNCT__
683 #define __FUNCT__ "PCBDDCSubSchursCreate"
684 PetscErrorCode PCBDDCSubSchursCreate(PCBDDCSubSchurs *sub_schurs)
685 {
686   PCBDDCSubSchurs schurs_ctx;
687   PetscErrorCode  ierr;
688 
689   PetscFunctionBegin;
690   ierr = PetscNew(&schurs_ctx);CHKERRQ(ierr);
691   schurs_ctx->n_subs = 0;
692   *sub_schurs = schurs_ctx;
693   PetscFunctionReturn(0);
694 }
695 
696 #undef __FUNCT__
697 #define __FUNCT__ "PCBDDCSubSchursDestroy"
698 PetscErrorCode PCBDDCSubSchursDestroy(PCBDDCSubSchurs *sub_schurs)
699 {
700   PetscErrorCode ierr;
701 
702   PetscFunctionBegin;
703   ierr = PCBDDCSubSchursReset(*sub_schurs);CHKERRQ(ierr);
704   ierr = PetscFree(*sub_schurs);CHKERRQ(ierr);
705   PetscFunctionReturn(0);
706 }
707 
708 #undef __FUNCT__
709 #define __FUNCT__ "PCBDDCSubSchursReset"
710 PetscErrorCode PCBDDCSubSchursReset(PCBDDCSubSchurs sub_schurs)
711 {
712   PetscInt       i;
713   PetscErrorCode ierr;
714 
715   PetscFunctionBegin;
716   ierr = MatDestroy(&sub_schurs->A);CHKERRQ(ierr);
717   ierr = MatDestroy(&sub_schurs->S);CHKERRQ(ierr);
718   ierr = ISDestroy(&sub_schurs->is_I);CHKERRQ(ierr);
719   ierr = ISDestroy(&sub_schurs->is_B);CHKERRQ(ierr);
720   ierr = ISLocalToGlobalMappingDestroy(&sub_schurs->l2gmap);CHKERRQ(ierr);
721   ierr = ISLocalToGlobalMappingDestroy(&sub_schurs->BtoNmap);CHKERRQ(ierr);
722   ierr = MatDestroy(&sub_schurs->S_Ej_all);CHKERRQ(ierr);
723   ierr = MatDestroy(&sub_schurs->sum_S_Ej_all);CHKERRQ(ierr);
724   ierr = MatDestroy(&sub_schurs->sum_S_Ej_tilda_all);CHKERRQ(ierr);
725   ierr = ISDestroy(&sub_schurs->is_Ej_all);CHKERRQ(ierr);
726   ierr = ISDestroy(&sub_schurs->is_Ej_com);CHKERRQ(ierr);
727   for (i=0;i<sub_schurs->n_subs;i++) {
728     ierr = ISDestroy(&sub_schurs->is_subs[i]);CHKERRQ(ierr);
729     ierr = MatDestroy(&sub_schurs->S_Ej[i]);CHKERRQ(ierr);
730     ierr = MatDestroy(&sub_schurs->S_Ej_tilda[i]);CHKERRQ(ierr);
731   }
732   ierr = ISDestroy(&sub_schurs->is_I_layer);CHKERRQ(ierr);
733   if (sub_schurs->n_subs) {
734     ierr = PetscFree(sub_schurs->is_subs);CHKERRQ(ierr);
735     ierr = PetscFree2(sub_schurs->S_Ej,sub_schurs->S_Ej_tilda);CHKERRQ(ierr);
736     ierr = PetscFree2(sub_schurs->index_sequential,sub_schurs->index_parallel);CHKERRQ(ierr);
737     ierr = PetscFree(sub_schurs->auxglobal_sequential);CHKERRQ(ierr);
738     ierr = PetscFree(sub_schurs->auxglobal_parallel);CHKERRQ(ierr);
739   }
740   sub_schurs->n_subs = 0;
741   PetscFunctionReturn(0);
742 }
743 
744 #undef __FUNCT__
745 #define __FUNCT__ "PCBDDCAdjGetNextLayer_Private"
746 PETSC_STATIC_INLINE PetscErrorCode PCBDDCAdjGetNextLayer_Private(PetscInt* queue_tip,PetscInt n_prev,PetscBT touched,PetscInt* xadj,PetscInt* adjncy,PetscInt* n_added)
747 {
748   PetscInt       i,j,n;
749   PetscErrorCode ierr;
750 
751   PetscFunctionBegin;
752   n = 0;
753   for (i=-n_prev;i<0;i++) {
754     PetscInt start_dof = queue_tip[i];
755     for (j=xadj[start_dof];j<xadj[start_dof+1];j++) {
756       PetscInt dof = adjncy[j];
757       if (!PetscBTLookup(touched,dof)) {
758         ierr = PetscBTSet(touched,dof);CHKERRQ(ierr);
759         queue_tip[n] = dof;
760         n++;
761       }
762     }
763   }
764   *n_added = n;
765   PetscFunctionReturn(0);
766 }
767