xref: /petsc/src/ksp/pc/impls/bddc/bddcschurs.c (revision b96c3477251e15db97d753dc3617251fb54656ca) 
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)
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   PetscBool              compute_Stilda=PETSC_FALSE;
101   PetscErrorCode         ierr;
102 
103   PetscFunctionBegin;
104   /* get Schur complement matrices */
105   if (!sub_schurs->use_mumps) {
106     ierr = MatSchurComplementGetSubMatrices(sub_schurs->S,&A_II,NULL,&A_IB,&A_BI,&A_BB);CHKERRQ(ierr);
107     ierr = PetscMalloc4(sub_schurs->n_subs,&is_subset_B,
108                         sub_schurs->n_subs,&AE_IE,
109                         sub_schurs->n_subs,&AE_EI,
110                         sub_schurs->n_subs,&AE_EE);CHKERRQ(ierr);
111   }
112 
113   /* determine interior problems */
114   ierr = ISDestroy(&sub_schurs->is_I_layer);CHKERRQ(ierr);
115   if (nlayers >= 0 && xadj != NULL && adjncy != NULL) { /* Interior problems can be different from the original one */
116     PetscBT                touched;
117     const PetscInt*        idx_B;
118     PetscInt               n_I,n_B,n_local_dofs,n_prev_added,j,layer,*local_numbering;
119 
120     /* get sizes */
121     ierr = ISGetLocalSize(sub_schurs->is_I,&n_I);CHKERRQ(ierr);
122     ierr = ISGetLocalSize(sub_schurs->is_B,&n_B);CHKERRQ(ierr);
123 
124     ierr = PetscMalloc1(n_I+n_B,&local_numbering);CHKERRQ(ierr);
125     ierr = PetscBTCreate(n_I+n_B,&touched);CHKERRQ(ierr);
126     ierr = PetscBTMemzero(n_I+n_B,touched);CHKERRQ(ierr);
127 
128     /* all boundary dofs must be skipped when adding layers */
129     ierr = ISGetIndices(sub_schurs->is_B,&idx_B);CHKERRQ(ierr);
130     for (j=0;j<n_B;j++) {
131       ierr = PetscBTSet(touched,idx_B[j]);CHKERRQ(ierr);
132     }
133     ierr = PetscMemcpy(local_numbering,idx_B,n_B*sizeof(PetscInt));CHKERRQ(ierr);
134     ierr = ISRestoreIndices(sub_schurs->is_B,&idx_B);CHKERRQ(ierr);
135 
136     /* add prescribed number of layers of dofs */
137     n_local_dofs = n_B;
138     n_prev_added = n_B;
139     for (layer=0;layer<nlayers;layer++) {
140       PetscInt n_added;
141       if (n_local_dofs == n_I+n_B) break;
142       if (n_local_dofs > n_I+n_B) {
143         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);
144       }
145       ierr = PCBDDCAdjGetNextLayer_Private(local_numbering+n_local_dofs,n_prev_added,touched,xadj,adjncy,&n_added);CHKERRQ(ierr);
146       n_prev_added = n_added;
147       n_local_dofs += n_added;
148       if (!n_added) break;
149     }
150     ierr = PetscBTDestroy(&touched);CHKERRQ(ierr);
151 
152     /* IS for I layer dofs in original numbering */
153     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);
154     ierr = PetscFree(local_numbering);CHKERRQ(ierr);
155     ierr = ISSort(sub_schurs->is_I_layer);CHKERRQ(ierr);
156     /* IS for I layer dofs in I numbering */
157     if (!sub_schurs->use_mumps) {
158       ISLocalToGlobalMapping ItoNmap;
159       ierr = ISLocalToGlobalMappingCreateIS(sub_schurs->is_I,&ItoNmap);CHKERRQ(ierr);
160       ierr = ISGlobalToLocalMappingApplyIS(ItoNmap,IS_GTOLM_DROP,sub_schurs->is_I_layer,&is_I);CHKERRQ(ierr);
161       ierr = ISLocalToGlobalMappingDestroy(&ItoNmap);CHKERRQ(ierr);
162 
163       /* II block */
164       ierr = MatGetSubMatrix(A_II,is_I,is_I,MAT_INITIAL_MATRIX,&AE_II);CHKERRQ(ierr);
165     }
166   } else {
167     PetscInt n_I;
168 
169     /* IS for I dofs in original numbering */
170     ierr = PetscObjectReference((PetscObject)sub_schurs->is_I);CHKERRQ(ierr);
171     sub_schurs->is_I_layer = sub_schurs->is_I;
172 
173     /* IS for I dofs in I numbering (strided 1) */
174     if (!sub_schurs->use_mumps) {
175       ierr = ISGetSize(sub_schurs->is_I,&n_I);CHKERRQ(ierr);
176       ierr = ISCreateStride(PetscObjectComm((PetscObject)sub_schurs->is_I),n_I,0,1,&is_I);CHKERRQ(ierr);
177 
178       /* II block is the same */
179       ierr = PetscObjectReference((PetscObject)A_II);CHKERRQ(ierr);
180       AE_II = A_II;
181     }
182   }
183 
184   /* Get info on subset sizes and sum of all subsets sizes */
185   max_subset_size = 0;
186   local_size = 0;
187   for (i=0;i<sub_schurs->n_subs_seq;i++) {
188     PetscInt j = sub_schurs->index_sequential[i];
189     ierr = ISGetLocalSize(sub_schurs->is_subs[j],&subset_size);CHKERRQ(ierr);
190     max_subset_size = PetscMax(subset_size,max_subset_size);
191     local_size += subset_size;
192   }
193 
194   /* Work arrays for local indices */
195   ierr = PetscMalloc1(local_size,&all_local_idx_B);CHKERRQ(ierr);
196   extra = 0;
197   if (sub_schurs->use_mumps) {
198     ierr = ISGetLocalSize(sub_schurs->is_I_layer,&extra);CHKERRQ(ierr);
199   }
200   ierr = PetscMalloc1(local_size+extra,&all_local_idx_N);CHKERRQ(ierr);
201   if (extra) {
202     const PetscInt *idxs;
203     ierr = ISGetIndices(sub_schurs->is_I_layer,&idxs);CHKERRQ(ierr);
204     ierr = PetscMemcpy(all_local_idx_N,idxs,extra*sizeof(PetscInt));CHKERRQ(ierr);
205     ierr = ISRestoreIndices(sub_schurs->is_I_layer,&idxs);CHKERRQ(ierr);
206   }
207   ierr = PetscMalloc1(local_size,&nnz);CHKERRQ(ierr);
208 
209   /* Get local indices in local numbering */
210   local_size = 0;
211   for (i=0;i<sub_schurs->n_subs_seq;i++) {
212     PetscInt j;
213     const    PetscInt *idxs;
214 
215     PetscInt local_problem_index = sub_schurs->index_sequential[i];
216     ierr = ISGetLocalSize(sub_schurs->is_subs[local_problem_index],&subset_size);CHKERRQ(ierr);
217     ierr = ISGetIndices(sub_schurs->is_subs[local_problem_index],&idxs);CHKERRQ(ierr);
218     /* subset indices in local numbering */
219     ierr = PetscMemcpy(all_local_idx_N+local_size+extra,idxs,subset_size*sizeof(PetscInt));CHKERRQ(ierr);
220     ierr = ISRestoreIndices(sub_schurs->is_subs[local_problem_index],&idxs);CHKERRQ(ierr);
221     for (j=0;j<subset_size;j++) nnz[local_size+j] = subset_size;
222     local_size += subset_size;
223   }
224 
225   S_all = NULL;
226   if (!sub_schurs->use_mumps) {
227     /* subsets in original and boundary numbering */
228     for (i=0;i<sub_schurs->n_subs;i++) {
229       ierr = ISGlobalToLocalMappingApplyIS(sub_schurs->BtoNmap,IS_GTOLM_DROP,sub_schurs->is_subs[i],&is_subset_B[i]);CHKERRQ(ierr);
230     }
231 
232     /* EE block */
233     for (i=0;i<sub_schurs->n_subs;i++) {
234       ierr = MatGetSubMatrix(A_BB,is_subset_B[i],is_subset_B[i],MAT_INITIAL_MATRIX,&AE_EE[i]);CHKERRQ(ierr);
235     }
236     /* IE block */
237     for (i=0;i<sub_schurs->n_subs;i++) {
238       ierr = MatGetSubMatrix(A_IB,is_I,is_subset_B[i],MAT_INITIAL_MATRIX,&AE_IE[i]);CHKERRQ(ierr);
239     }
240     /* EI block */
241     for (i=0;i<sub_schurs->n_subs;i++) {
242       ierr = MatGetSubMatrix(A_BI,is_subset_B[i],is_I,MAT_INITIAL_MATRIX,&AE_EI[i]);CHKERRQ(ierr);
243     }
244 
245     /* setup Schur complements on subset */
246     for (i=0;i<sub_schurs->n_subs;i++) {
247       ierr = MatDestroy(&sub_schurs->S_Ej[i]);CHKERRQ(ierr);
248       ierr = MatCreateSchurComplement(AE_II,AE_II,AE_IE[i],AE_EI[i],AE_EE[i],&sub_schurs->S_Ej[i]);CHKERRQ(ierr);
249       if (AE_II == A_II) { /* we can reuse the same ksp */
250         KSP ksp;
251         ierr = MatSchurComplementGetKSP(sub_schurs->S,&ksp);CHKERRQ(ierr);
252         ierr = MatSchurComplementSetKSP(sub_schurs->S_Ej[i],ksp);CHKERRQ(ierr);
253       } else { /* build new ksp object which inherits ksp and pc types from the original one */
254         KSP      origksp,schurksp;
255         PC       origpc,schurpc;
256         KSPType  ksp_type;
257         PCType   pc_type;
258         PetscInt n_internal;
259 
260         ierr = MatSchurComplementGetKSP(sub_schurs->S,&origksp);CHKERRQ(ierr);
261         ierr = MatSchurComplementGetKSP(sub_schurs->S_Ej[i],&schurksp);CHKERRQ(ierr);
262         ierr = KSPGetType(origksp,&ksp_type);CHKERRQ(ierr);
263         ierr = KSPSetType(schurksp,ksp_type);CHKERRQ(ierr);
264         ierr = KSPGetPC(schurksp,&schurpc);CHKERRQ(ierr);
265         ierr = KSPGetPC(origksp,&origpc);CHKERRQ(ierr);
266         ierr = PCGetType(origpc,&pc_type);CHKERRQ(ierr);
267         ierr = PCSetType(schurpc,pc_type);CHKERRQ(ierr);
268         ierr = ISGetSize(is_I,&n_internal);CHKERRQ(ierr);
269         if (n_internal) { /* UMFPACK gives error with 0 sized problems */
270           MatSolverPackage solver=NULL;
271           ierr = PCFactorGetMatSolverPackage(origpc,(const MatSolverPackage*)&solver);CHKERRQ(ierr);
272           if (solver) {
273             ierr = PCFactorSetMatSolverPackage(schurpc,solver);CHKERRQ(ierr);
274           }
275         }
276         ierr = KSPSetUp(schurksp);CHKERRQ(ierr);
277       }
278     }
279     /* free */
280     ierr = ISDestroy(&is_I);CHKERRQ(ierr);
281     ierr = MatDestroy(&AE_II);CHKERRQ(ierr);
282     for (i=0;i<sub_schurs->n_subs;i++) {
283       ierr = MatDestroy(&AE_EE[i]);CHKERRQ(ierr);
284       ierr = MatDestroy(&AE_IE[i]);CHKERRQ(ierr);
285       ierr = MatDestroy(&AE_EI[i]);CHKERRQ(ierr);
286       ierr = ISDestroy(&is_subset_B[i]);CHKERRQ(ierr);
287     }
288     ierr = PetscFree4(is_subset_B,AE_IE,AE_EI,AE_EE);CHKERRQ(ierr);
289 #if defined(PETSC_HAVE_MUMPS)
290   } else {
291     Mat           A,F;
292     IS            is_A_all;
293     PetscBool     is_symmetric;
294     PetscInt      *idxs_schur,n_I;
295 
296     /* get working mat */
297     ierr = ISGetLocalSize(sub_schurs->is_I_layer,&n_I);CHKERRQ(ierr);
298     ierr = ISCreateGeneral(PETSC_COMM_SELF,local_size+n_I,all_local_idx_N,PETSC_COPY_VALUES,&is_A_all);CHKERRQ(ierr);
299     ierr = MatGetSubMatrixUnsorted(sub_schurs->A,is_A_all,is_A_all,MAT_INITIAL_MATRIX,&A);CHKERRQ(ierr);
300     ierr = ISDestroy(&is_A_all);CHKERRQ(ierr);
301 
302     ierr = MatIsSymmetric(sub_schurs->A,0.0,&is_symmetric);CHKERRQ(ierr);
303     if (is_symmetric) {
304       ierr = MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_CHOLESKY,&F);CHKERRQ(ierr);
305     } else {
306       ierr = MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_LU,&F);CHKERRQ(ierr);
307     }
308 
309     /* subsets ordered last */
310     ierr = PetscMalloc1(local_size,&idxs_schur);CHKERRQ(ierr);
311     for (i=0;i<local_size;i++) {
312       idxs_schur[i] = n_I+i+1;
313     }
314     ierr = MatMumpsSetSchurIndices(F,local_size,idxs_schur);CHKERRQ(ierr);
315     ierr = PetscFree(idxs_schur);CHKERRQ(ierr);
316 
317     /* factorization step */
318     if (is_symmetric) {
319       ierr = MatCholeskyFactorSymbolic(F,A,NULL,NULL);CHKERRQ(ierr);
320       ierr = MatCholeskyFactorNumeric(F,A,NULL);CHKERRQ(ierr);
321     } else {
322       ierr = MatLUFactorSymbolic(F,A,NULL,NULL,NULL);CHKERRQ(ierr);
323       ierr = MatLUFactorNumeric(F,A,NULL);CHKERRQ(ierr);
324     }
325 
326     /* get explicit Schur Complement computed during numeric factorization */
327     ierr = MatMumpsGetSchurComplement(F,&S_all);CHKERRQ(ierr);
328 
329     /* free workspace */
330     ierr = MatDestroy(&A);CHKERRQ(ierr);
331     ierr = MatDestroy(&F);CHKERRQ(ierr);
332 
333 #endif
334   }
335 
336   if (!sub_schurs->n_subs_seq_g) {
337     PetscFunctionReturn(0);
338   }
339 
340   /* subset indices in local boundary numbering */
341   if (!sub_schurs->is_Ej_all) {
342     ierr = ISGlobalToLocalMappingApply(sub_schurs->BtoNmap,IS_GTOLM_DROP,local_size,all_local_idx_N+extra,&subset_size,all_local_idx_B);CHKERRQ(ierr);
343     if (subset_size != local_size) {
344       SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in sub_schurs serial (BtoNmap)! %d != %d\n",subset_size,local_size);
345     }
346     ierr = ISCreateGeneral(PETSC_COMM_SELF,local_size,all_local_idx_B,PETSC_OWN_POINTER,&sub_schurs->is_Ej_all);CHKERRQ(ierr);
347   }
348 
349   /* Local matrix of all local Schur on subsets transposed */
350   if (!sub_schurs->S_Ej_all) {
351     ierr = MatCreate(PETSC_COMM_SELF,&sub_schurs->S_Ej_all);CHKERRQ(ierr);
352     ierr = MatSetSizes(sub_schurs->S_Ej_all,PETSC_DECIDE,PETSC_DECIDE,local_size,local_size);CHKERRQ(ierr);
353     ierr = MatSetType(sub_schurs->S_Ej_all,MATAIJ);CHKERRQ(ierr);
354     ierr = MatSeqAIJSetPreallocation(sub_schurs->S_Ej_all,0,nnz);CHKERRQ(ierr);
355   } else {
356     ierr = MatZeroEntries(sub_schurs->S_Ej_all);CHKERRQ(ierr);
357   }
358 
359   S_Ej_tilda_all = 0;
360   if (!sub_schurs->use_mumps) {
361     PetscScalar *fill_vals;
362     PetscInt    *dummy_idx;
363 
364     /* Work arrays */
365     ierr = PetscMalloc1(max_subset_size,&dummy_idx);CHKERRQ(ierr);
366 
367     /* Loop on local problems end compute Schur complements explicitly */
368     local_size = 0;
369     for (i=0;i<sub_schurs->n_subs_seq;i++) {
370       Mat       S_Ej_expl;
371       PetscInt  j,lpi = sub_schurs->index_sequential[i];
372       PetscBool Sdense;
373 
374       ierr = PCBDDCComputeExplicitSchur(sub_schurs->S_Ej[lpi],&S_Ej_expl);CHKERRQ(ierr);
375       ierr = ISGetLocalSize(sub_schurs->is_subs[lpi],&subset_size);CHKERRQ(ierr);
376       ierr = PetscObjectTypeCompare((PetscObject)S_Ej_expl,MATSEQDENSE,&Sdense);CHKERRQ(ierr);
377       if (Sdense) {
378         for (j=0;j<subset_size;j++) {
379           dummy_idx[j]=local_size+j;
380         }
381         ierr = MatDenseGetArray(S_Ej_expl,&fill_vals);CHKERRQ(ierr);
382         ierr = MatSetValues(sub_schurs->S_Ej_all,subset_size,dummy_idx,subset_size,dummy_idx,fill_vals,INSERT_VALUES);CHKERRQ(ierr);
383         ierr = MatDenseRestoreArray(S_Ej_expl,&fill_vals);CHKERRQ(ierr);
384       } else {
385         SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not yet implemented for sparse matrices");
386       }
387       local_size += subset_size;
388       ierr = MatDestroy(&sub_schurs->S_Ej[lpi]);CHKERRQ(ierr);
389       sub_schurs->S_Ej[lpi] = S_Ej_expl;
390     }
391     /* Stildas are not computed without mumps */
392     for (i=0;i<sub_schurs->n_subs;i++) {
393       ierr = MatDestroy(&sub_schurs->S_Ej_tilda[i]);CHKERRQ(ierr);
394     }
395     ierr = PetscFree(dummy_idx);CHKERRQ(ierr);
396   } else {
397     PetscInt  *dummy_idx,n_all;
398     PetscBool is_symmetric;
399 
400     compute_Stilda=PETSC_TRUE;
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;
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,NULL);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   for (i=0;i<sub_schurs->n_subs;i++) {
656     ierr = ISSort(sub_schurs->is_subs[i]);CHKERRQ(ierr);
657   }
658   sub_schurs->index_sequential = index_sequential;
659   sub_schurs->index_parallel = index_parallel;
660   sub_schurs->auxglobal_sequential = auxglobal_sequential;
661   sub_schurs->auxglobal_parallel = auxglobal_parallel;
662 
663 
664   /* allocate space for schur complements */
665   ierr = PetscCalloc2(sub_schurs->n_subs,&sub_schurs->S_Ej,
666                       sub_schurs->n_subs,&sub_schurs->S_Ej_tilda);CHKERRQ(ierr);
667   sub_schurs->S_Ej_all = NULL;
668   sub_schurs->sum_S_Ej_all = NULL;
669   sub_schurs->sum_S_Ej_tilda_all = NULL;
670   sub_schurs->is_Ej_all = NULL;
671 
672   /* free workspace */
673   ierr = PetscFree(auxmapping);CHKERRQ(ierr);
674   ierr = PetscFree2(auxlocal_sequential,auxlocal_parallel);CHKERRQ(ierr);
675 
676   PetscFunctionReturn(0);
677 }
678 
679 #undef __FUNCT__
680 #define __FUNCT__ "PCBDDCSubSchursCreate"
681 PetscErrorCode PCBDDCSubSchursCreate(PCBDDCSubSchurs *sub_schurs)
682 {
683   PCBDDCSubSchurs schurs_ctx;
684   PetscErrorCode  ierr;
685 
686   PetscFunctionBegin;
687   ierr = PetscNew(&schurs_ctx);CHKERRQ(ierr);
688   schurs_ctx->n_subs = 0;
689   *sub_schurs = schurs_ctx;
690   PetscFunctionReturn(0);
691 }
692 
693 #undef __FUNCT__
694 #define __FUNCT__ "PCBDDCSubSchursDestroy"
695 PetscErrorCode PCBDDCSubSchursDestroy(PCBDDCSubSchurs *sub_schurs)
696 {
697   PetscErrorCode ierr;
698 
699   PetscFunctionBegin;
700   ierr = PCBDDCSubSchursReset(*sub_schurs);CHKERRQ(ierr);
701   ierr = PetscFree(*sub_schurs);CHKERRQ(ierr);
702   PetscFunctionReturn(0);
703 }
704 
705 #undef __FUNCT__
706 #define __FUNCT__ "PCBDDCSubSchursReset"
707 PetscErrorCode PCBDDCSubSchursReset(PCBDDCSubSchurs sub_schurs)
708 {
709   PetscInt       i;
710   PetscErrorCode ierr;
711 
712   PetscFunctionBegin;
713   ierr = MatDestroy(&sub_schurs->A);CHKERRQ(ierr);
714   ierr = MatDestroy(&sub_schurs->S);CHKERRQ(ierr);
715   ierr = ISDestroy(&sub_schurs->is_I);CHKERRQ(ierr);
716   ierr = ISDestroy(&sub_schurs->is_B);CHKERRQ(ierr);
717   ierr = ISLocalToGlobalMappingDestroy(&sub_schurs->l2gmap);CHKERRQ(ierr);
718   ierr = ISLocalToGlobalMappingDestroy(&sub_schurs->BtoNmap);CHKERRQ(ierr);
719   ierr = MatDestroy(&sub_schurs->S_Ej_all);CHKERRQ(ierr);
720   ierr = MatDestroy(&sub_schurs->sum_S_Ej_all);CHKERRQ(ierr);
721   ierr = MatDestroy(&sub_schurs->sum_S_Ej_tilda_all);CHKERRQ(ierr);
722   ierr = ISDestroy(&sub_schurs->is_Ej_all);CHKERRQ(ierr);
723   for (i=0;i<sub_schurs->n_subs;i++) {
724     ierr = ISDestroy(&sub_schurs->is_subs[i]);CHKERRQ(ierr);
725     ierr = MatDestroy(&sub_schurs->S_Ej[i]);CHKERRQ(ierr);
726     ierr = MatDestroy(&sub_schurs->S_Ej_tilda[i]);CHKERRQ(ierr);
727   }
728   ierr = ISDestroy(&sub_schurs->is_I_layer);CHKERRQ(ierr);
729   if (sub_schurs->n_subs) {
730     ierr = PetscFree(sub_schurs->is_subs);CHKERRQ(ierr);
731     ierr = PetscFree2(sub_schurs->S_Ej,sub_schurs->S_Ej_tilda);CHKERRQ(ierr);
732     ierr = PetscFree2(sub_schurs->index_sequential,sub_schurs->index_parallel);CHKERRQ(ierr);
733     ierr = PetscFree(sub_schurs->auxglobal_sequential);CHKERRQ(ierr);
734     ierr = PetscFree(sub_schurs->auxglobal_parallel);CHKERRQ(ierr);
735   }
736   sub_schurs->n_subs = 0;
737   PetscFunctionReturn(0);
738 }
739 
740 #undef __FUNCT__
741 #define __FUNCT__ "PCBDDCAdjGetNextLayer_Private"
742 PETSC_STATIC_INLINE PetscErrorCode PCBDDCAdjGetNextLayer_Private(PetscInt* queue_tip,PetscInt n_prev,PetscBT touched,PetscInt* xadj,PetscInt* adjncy,PetscInt* n_added)
743 {
744   PetscInt       i,j,n;
745   PetscErrorCode ierr;
746 
747   PetscFunctionBegin;
748   n = 0;
749   for (i=-n_prev;i<0;i++) {
750     PetscInt start_dof = queue_tip[i];
751     for (j=xadj[start_dof];j<xadj[start_dof+1];j++) {
752       PetscInt dof = adjncy[j];
753       if (!PetscBTLookup(touched,dof)) {
754         ierr = PetscBTSet(touched,dof);CHKERRQ(ierr);
755         queue_tip[n] = dof;
756         n++;
757       }
758     }
759   }
760   *n_added = n;
761   PetscFunctionReturn(0);
762 }
763