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