xref: /petsc/src/ksp/pc/impls/bddc/bddcscalingbasic.c (revision 95dccacacae8a8fc0b691f9b4fba69a249b61188) 
1 #include <../src/ksp/pc/impls/bddc/bddc.h>
2 #include <../src/ksp/pc/impls/bddc/bddcprivate.h>
3 
4 /* prototypes for deluxe functions */
5 static PetscErrorCode PCBDDCScalingCreate_Deluxe(PC);
6 static PetscErrorCode PCBDDCScalingDestroy_Deluxe(PC);
7 static PetscErrorCode PCBDDCScalingSetUp_Deluxe(PC);
8 static PetscErrorCode PCBDDCScalingSetUp_Deluxe_Private(PC);
9 static PetscErrorCode PCBDDCScalingReset_Deluxe_Solvers(PCBDDCDeluxeScaling);
10 
11 #undef __FUNCT__
12 #define __FUNCT__ "PCBDDCScalingExtension_Basic"
13 static PetscErrorCode PCBDDCScalingExtension_Basic(PC pc, Vec local_interface_vector, Vec global_vector)
14 {
15   PC_IS*         pcis = (PC_IS*)pc->data;
16   PC_BDDC*       pcbddc = (PC_BDDC*)pc->data;
17   PetscErrorCode ierr;
18 
19   PetscFunctionBegin;
20   /* Apply partition of unity */
21   ierr = VecPointwiseMult(pcbddc->work_scaling,pcis->D,local_interface_vector);CHKERRQ(ierr);
22   ierr = VecSet(global_vector,0.0);CHKERRQ(ierr);
23   ierr = VecScatterBegin(pcis->global_to_B,pcbddc->work_scaling,global_vector,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
24   ierr = VecScatterEnd(pcis->global_to_B,pcbddc->work_scaling,global_vector,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
25   PetscFunctionReturn(0);
26 }
27 
28 #undef __FUNCT__
29 #define __FUNCT__ "PCBDDCScalingExtension_Deluxe"
30 static PetscErrorCode PCBDDCScalingExtension_Deluxe(PC pc, Vec x, Vec y)
31 {
32   PC_IS*              pcis=(PC_IS*)pc->data;
33   PC_BDDC*            pcbddc=(PC_BDDC*)pc->data;
34   PCBDDCDeluxeScaling deluxe_ctx = pcbddc->deluxe_ctx;
35   PetscErrorCode      ierr;
36 
37   PetscFunctionBegin;
38   ierr = VecSet(pcbddc->work_scaling,0.0);CHKERRQ(ierr);
39   ierr = VecSet(y,0.0);CHKERRQ(ierr);
40   if (deluxe_ctx->n_simple) { /* scale deluxe vertices using diagonal scaling */
41     PetscInt          i;
42     const PetscScalar *array_x,*array_D;
43     PetscScalar       *array;
44     ierr = VecGetArrayRead(x,&array_x);CHKERRQ(ierr);
45     ierr = VecGetArrayRead(pcis->D,&array_D);CHKERRQ(ierr);
46     ierr = VecGetArray(pcbddc->work_scaling,&array);CHKERRQ(ierr);
47     for (i=0;i<deluxe_ctx->n_simple;i++) {
48       array[deluxe_ctx->idx_simple_B[i]] = array_x[deluxe_ctx->idx_simple_B[i]]*array_D[deluxe_ctx->idx_simple_B[i]];
49     }
50     ierr = VecRestoreArray(pcbddc->work_scaling,&array);CHKERRQ(ierr);
51     ierr = VecRestoreArrayRead(pcis->D,&array_D);CHKERRQ(ierr);
52     ierr = VecRestoreArrayRead(x,&array_x);CHKERRQ(ierr);
53   }
54   /* sequential part : all problems and Schur applications collapsed into a single matrix vector multiplication or a matvec and a solve */
55   if (deluxe_ctx->seq_mat) {
56     ierr = VecScatterBegin(deluxe_ctx->seq_scctx,x,deluxe_ctx->seq_work1,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
57     ierr = VecScatterEnd(deluxe_ctx->seq_scctx,x,deluxe_ctx->seq_work1,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
58     ierr = MatMultTranspose(deluxe_ctx->seq_mat,deluxe_ctx->seq_work1,deluxe_ctx->seq_work2);CHKERRQ(ierr);
59     if (deluxe_ctx->seq_ksp) {
60       ierr = KSPSolveTranspose(deluxe_ctx->seq_ksp,deluxe_ctx->seq_work2,deluxe_ctx->seq_work2);CHKERRQ(ierr);
61     }
62     ierr = VecScatterBegin(deluxe_ctx->seq_scctx,deluxe_ctx->seq_work2,pcbddc->work_scaling,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
63     ierr = VecScatterEnd(deluxe_ctx->seq_scctx,deluxe_ctx->seq_work2,pcbddc->work_scaling,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
64   }
65   /* put local boundary part in global vector */
66   ierr = VecScatterBegin(pcis->global_to_B,pcbddc->work_scaling,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
67   ierr = VecScatterEnd(pcis->global_to_B,pcbddc->work_scaling,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
68   PetscFunctionReturn(0);
69 }
70 
71 #undef __FUNCT__
72 #define __FUNCT__ "PCBDDCScalingExtension"
73 PetscErrorCode PCBDDCScalingExtension(PC pc, Vec local_interface_vector, Vec global_vector)
74 {
75   PC_BDDC        *pcbddc=(PC_BDDC*)pc->data;
76   PetscErrorCode ierr;
77 
78   PetscFunctionBegin;
79   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
80   PetscValidHeaderSpecific(local_interface_vector,VEC_CLASSID,2);
81   PetscValidHeaderSpecific(global_vector,VEC_CLASSID,3);
82   if (local_interface_vector == pcbddc->work_scaling) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Local vector cannot be pcbddc->work_scaling!\n");
83   ierr = PetscUseMethod(pc,"PCBDDCScalingExtension_C",(PC,Vec,Vec),(pc,local_interface_vector,global_vector));CHKERRQ(ierr);
84   PetscFunctionReturn(0);
85 }
86 
87 #undef __FUNCT__
88 #define __FUNCT__ "PCBDDCScalingRestriction_Basic"
89 static PetscErrorCode PCBDDCScalingRestriction_Basic(PC pc, Vec global_vector, Vec local_interface_vector)
90 {
91   PetscErrorCode ierr;
92   PC_IS          *pcis = (PC_IS*)pc->data;
93 
94   PetscFunctionBegin;
95   ierr = VecScatterBegin(pcis->global_to_B,global_vector,local_interface_vector,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
96   ierr = VecScatterEnd(pcis->global_to_B,global_vector,local_interface_vector,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
97   /* Apply partition of unity */
98   ierr = VecPointwiseMult(local_interface_vector,pcis->D,local_interface_vector);CHKERRQ(ierr);
99   PetscFunctionReturn(0);
100 }
101 
102 #undef __FUNCT__
103 #define __FUNCT__ "PCBDDCScalingRestriction_Deluxe"
104 static PetscErrorCode PCBDDCScalingRestriction_Deluxe(PC pc, Vec x, Vec y)
105 {
106   PC_IS*              pcis=(PC_IS*)pc->data;
107   PC_BDDC*            pcbddc=(PC_BDDC*)pc->data;
108   PCBDDCDeluxeScaling deluxe_ctx = pcbddc->deluxe_ctx;
109   PetscErrorCode      ierr;
110 
111   PetscFunctionBegin;
112   /* get local boundary part of global vector */
113   ierr = VecScatterBegin(pcis->global_to_B,x,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
114   ierr = VecScatterEnd(pcis->global_to_B,x,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
115   if (deluxe_ctx->n_simple) { /* scale deluxe vertices using diagonal scaling */
116     PetscInt          i;
117     PetscScalar       *array_y;
118     const PetscScalar *array_D;
119     ierr = VecGetArray(y,&array_y);CHKERRQ(ierr);
120     ierr = VecGetArrayRead(pcis->D,&array_D);CHKERRQ(ierr);
121     for (i=0;i<deluxe_ctx->n_simple;i++) {
122       array_y[deluxe_ctx->idx_simple_B[i]] *= array_D[deluxe_ctx->idx_simple_B[i]];
123     }
124     ierr = VecRestoreArrayRead(pcis->D,&array_D);CHKERRQ(ierr);
125     ierr = VecRestoreArray(y,&array_y);CHKERRQ(ierr);
126   }
127   /* sequential part : all problems and Schur applications collapsed into a single matrix vector multiplication and ksp solution */
128   if (deluxe_ctx->seq_mat) {
129     ierr = VecScatterBegin(deluxe_ctx->seq_scctx,y,deluxe_ctx->seq_work1,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
130     ierr = VecScatterEnd(deluxe_ctx->seq_scctx,y,deluxe_ctx->seq_work1,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
131     if (deluxe_ctx->seq_ksp) {
132       ierr = KSPSolve(deluxe_ctx->seq_ksp,deluxe_ctx->seq_work1,deluxe_ctx->seq_work1);CHKERRQ(ierr);
133     }
134     ierr = MatMult(deluxe_ctx->seq_mat,deluxe_ctx->seq_work1,deluxe_ctx->seq_work2);CHKERRQ(ierr);
135     ierr = VecScatterBegin(deluxe_ctx->seq_scctx,deluxe_ctx->seq_work2,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
136     ierr = VecScatterEnd(deluxe_ctx->seq_scctx,deluxe_ctx->seq_work2,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
137   }
138   PetscFunctionReturn(0);
139 }
140 
141 #undef __FUNCT__
142 #define __FUNCT__ "PCBDDCScalingRestriction"
143 PetscErrorCode PCBDDCScalingRestriction(PC pc, Vec global_vector, Vec local_interface_vector)
144 {
145   PC_BDDC        *pcbddc=(PC_BDDC*)pc->data;
146   PetscErrorCode ierr;
147 
148   PetscFunctionBegin;
149   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
150   PetscValidHeaderSpecific(global_vector,VEC_CLASSID,2);
151   PetscValidHeaderSpecific(local_interface_vector,VEC_CLASSID,3);
152   if (local_interface_vector == pcbddc->work_scaling) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Local vector cannot be pcbddc->work_scaling!\n");
153   ierr = PetscUseMethod(pc,"PCBDDCScalingRestriction_C",(PC,Vec,Vec),(pc,global_vector,local_interface_vector));CHKERRQ(ierr);
154   PetscFunctionReturn(0);
155 }
156 
157 #undef __FUNCT__
158 #define __FUNCT__ "PCBDDCScalingSetUp"
159 PetscErrorCode PCBDDCScalingSetUp(PC pc)
160 {
161   PC_IS*         pcis=(PC_IS*)pc->data;
162   PC_BDDC*       pcbddc=(PC_BDDC*)pc->data;
163   PetscErrorCode ierr;
164 
165   PetscFunctionBegin;
166   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
167   /* create work vector for the operator */
168   ierr = VecDestroy(&pcbddc->work_scaling);CHKERRQ(ierr);
169   ierr = VecDuplicate(pcis->vec1_B,&pcbddc->work_scaling);CHKERRQ(ierr);
170   /* always rebuild pcis->D */
171   if (pcis->use_stiffness_scaling) {
172     ierr = MatGetDiagonal(pcbddc->local_mat,pcis->vec1_N);CHKERRQ(ierr);
173     ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
174     ierr = VecScatterEnd(pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
175   }
176   ierr = VecCopy(pcis->D,pcis->vec1_B);CHKERRQ(ierr);
177   ierr = VecSet(pcis->vec1_global,0.0);CHKERRQ(ierr);
178   ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
179   ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
180   ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
181   ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
182   ierr = VecPointwiseDivide(pcis->D,pcis->D,pcis->vec1_B);CHKERRQ(ierr);
183   /* now setup */
184   if (pcbddc->use_deluxe_scaling) {
185     if (!pcbddc->deluxe_ctx) {
186       ierr = PCBDDCScalingCreate_Deluxe(pc);CHKERRQ(ierr);
187     }
188     ierr = PCBDDCScalingSetUp_Deluxe(pc);CHKERRQ(ierr);
189     ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingRestriction_C",PCBDDCScalingRestriction_Deluxe);CHKERRQ(ierr);
190     ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingExtension_C",PCBDDCScalingExtension_Deluxe);CHKERRQ(ierr);
191   } else {
192     ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingRestriction_C",PCBDDCScalingRestriction_Basic);CHKERRQ(ierr);
193     ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingExtension_C",PCBDDCScalingExtension_Basic);CHKERRQ(ierr);
194   }
195 
196   /* test */
197   if (pcbddc->dbg_flag) {
198     Vec         vec2_global;
199     PetscViewer viewer=pcbddc->dbg_viewer;
200     PetscReal   error;
201 
202     /* extension -> from local to parallel */
203     ierr = VecSet(pcis->vec1_global,0.0);CHKERRQ(ierr);
204     ierr = VecSetRandom(pcis->vec1_B,NULL);CHKERRQ(ierr);
205     ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
206     ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
207     ierr = VecDuplicate(pcis->vec1_global,&vec2_global);CHKERRQ(ierr);
208     ierr = VecCopy(pcis->vec1_global,vec2_global);CHKERRQ(ierr);
209 
210     ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
211     ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
212     ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,pcis->vec1_global);CHKERRQ(ierr);
213     ierr = VecAXPY(pcis->vec1_global,-1.0,vec2_global);CHKERRQ(ierr);
214     ierr = VecNorm(pcis->vec1_global,NORM_INFINITY,&error);CHKERRQ(ierr);
215     ierr = PetscViewerASCIIPrintf(viewer,"Error scaling extension %1.14e\n",error);CHKERRQ(ierr);
216     if (error>1.e-8 && pcbddc->dbg_flag>1) {
217       ierr = VecView(pcis->vec1_global,viewer);CHKERRQ(ierr);
218     }
219     ierr = VecDestroy(&vec2_global);CHKERRQ(ierr);
220 
221     /* restriction -> from parallel to local */
222     ierr = VecSet(pcis->vec1_global,0.0);CHKERRQ(ierr);
223     ierr = VecSetRandom(pcis->vec1_B,NULL);CHKERRQ(ierr);
224     ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
225     ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
226 
227     ierr = PCBDDCScalingRestriction(pc,pcis->vec1_global,pcis->vec1_B);CHKERRQ(ierr);
228     ierr = VecScale(pcis->vec1_B,-1.0);CHKERRQ(ierr);
229     ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
230     ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
231     ierr = VecNorm(pcis->vec1_global,NORM_INFINITY,&error);CHKERRQ(ierr);
232     ierr = PetscViewerASCIIPrintf(viewer,"Error scaling restriction %1.14e\n",error);CHKERRQ(ierr);
233     if (error>1.e-8 && pcbddc->dbg_flag>1) {
234       ierr = VecView(pcis->vec1_global,viewer);CHKERRQ(ierr);
235     }
236   }
237   PetscFunctionReturn(0);
238 }
239 
240 #undef __FUNCT__
241 #define __FUNCT__ "PCBDDCScalingDestroy"
242 PetscErrorCode PCBDDCScalingDestroy(PC pc)
243 {
244   PC_BDDC* pcbddc=(PC_BDDC*)pc->data;
245   PetscErrorCode ierr;
246 
247   PetscFunctionBegin;
248   if (pcbddc->deluxe_ctx) {
249     ierr = PCBDDCScalingDestroy_Deluxe(pc);CHKERRQ(ierr);
250   }
251   ierr = VecDestroy(&pcbddc->work_scaling);CHKERRQ(ierr);
252   /* remove functions */
253   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingRestriction_C",NULL);CHKERRQ(ierr);
254   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingExtension_C",NULL);CHKERRQ(ierr);
255   PetscFunctionReturn(0);
256 }
257 
258 #undef __FUNCT__
259 #define __FUNCT__ "PCBDDCScalingCreate_Deluxe"
260 static PetscErrorCode PCBDDCScalingCreate_Deluxe(PC pc)
261 {
262   PC_BDDC*            pcbddc=(PC_BDDC*)pc->data;
263   PCBDDCDeluxeScaling deluxe_ctx;
264   PetscErrorCode      ierr;
265 
266   PetscFunctionBegin;
267   ierr = PetscNew(&deluxe_ctx);CHKERRQ(ierr);
268   pcbddc->deluxe_ctx = deluxe_ctx;
269   PetscFunctionReturn(0);
270 }
271 
272 #undef __FUNCT__
273 #define __FUNCT__ "PCBDDCScalingDestroy_Deluxe"
274 static PetscErrorCode PCBDDCScalingDestroy_Deluxe(PC pc)
275 {
276   PC_BDDC*            pcbddc=(PC_BDDC*)pc->data;
277   PetscErrorCode      ierr;
278 
279   PetscFunctionBegin;
280   ierr = PCBDDCScalingReset_Deluxe_Solvers(pcbddc->deluxe_ctx);CHKERRQ(ierr);
281   ierr = PetscFree(pcbddc->deluxe_ctx);CHKERRQ(ierr);
282   PetscFunctionReturn(0);
283 }
284 
285 #undef __FUNCT__
286 #define __FUNCT__ "PCBDDCScalingReset_Deluxe_Solvers"
287 static PetscErrorCode PCBDDCScalingReset_Deluxe_Solvers(PCBDDCDeluxeScaling deluxe_ctx)
288 {
289   PetscErrorCode      ierr;
290 
291   PetscFunctionBegin;
292   ierr = PetscFree(deluxe_ctx->idx_simple_B);CHKERRQ(ierr);
293   deluxe_ctx->n_simple = 0;
294   ierr = VecScatterDestroy(&deluxe_ctx->seq_scctx);CHKERRQ(ierr);
295   ierr = VecDestroy(&deluxe_ctx->seq_work1);CHKERRQ(ierr);
296   ierr = VecDestroy(&deluxe_ctx->seq_work2);CHKERRQ(ierr);
297   ierr = MatDestroy(&deluxe_ctx->seq_mat);CHKERRQ(ierr);
298   ierr = KSPDestroy(&deluxe_ctx->seq_ksp);CHKERRQ(ierr);
299   PetscFunctionReturn(0);
300 }
301 
302 #undef __FUNCT__
303 #define __FUNCT__ "PCBDDCScalingSetUp_Deluxe"
304 static PetscErrorCode PCBDDCScalingSetUp_Deluxe(PC pc)
305 {
306   PC_IS               *pcis=(PC_IS*)pc->data;
307   PC_BDDC             *pcbddc=(PC_BDDC*)pc->data;
308   PCBDDCDeluxeScaling deluxe_ctx=pcbddc->deluxe_ctx;
309   PCBDDCSubSchurs     sub_schurs=pcbddc->sub_schurs;
310   PetscErrorCode      ierr;
311 
312   PetscFunctionBegin;
313   /* (TODO: reuse) throw away the solvers */
314   ierr = PCBDDCScalingReset_Deluxe_Solvers(deluxe_ctx);CHKERRQ(ierr);
315 
316   /* Compute data structures to solve sequential problems */
317   ierr = PCBDDCScalingSetUp_Deluxe_Private(pc);CHKERRQ(ierr);
318 
319   /* diagonal scaling on interface dofs not contained in cc */
320   if (sub_schurs->is_vertices || sub_schurs->is_dir) {
321     PetscInt n_com,n_dir;
322     n_com = 0;
323     if (sub_schurs->is_vertices) {
324       ierr = ISGetLocalSize(sub_schurs->is_vertices,&n_com);CHKERRQ(ierr);
325     }
326     n_dir = 0;
327     if (sub_schurs->is_dir) {
328       ierr = ISGetLocalSize(sub_schurs->is_dir,&n_dir);CHKERRQ(ierr);
329     }
330     deluxe_ctx->n_simple = n_dir + n_com;
331     ierr = PetscMalloc1(deluxe_ctx->n_simple,&deluxe_ctx->idx_simple_B);CHKERRQ(ierr);
332     if (sub_schurs->is_vertices) {
333       PetscInt       nmap;
334       const PetscInt *idxs;
335 
336       ierr = ISGetIndices(sub_schurs->is_vertices,&idxs);CHKERRQ(ierr);
337       ierr = ISGlobalToLocalMappingApply(pcis->BtoNmap,IS_GTOLM_DROP,n_com,idxs,&nmap,deluxe_ctx->idx_simple_B);CHKERRQ(ierr);
338       if (nmap != n_com) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error when mapping simply scaled dofs (is_vertices)! %D != %D",nmap,n_com);
339       ierr = ISRestoreIndices(sub_schurs->is_vertices,&idxs);CHKERRQ(ierr);
340     }
341     if (sub_schurs->is_dir) {
342       PetscInt       nmap;
343       const PetscInt *idxs;
344 
345       ierr = ISGetIndices(sub_schurs->is_dir,&idxs);CHKERRQ(ierr);
346       ierr = ISGlobalToLocalMappingApply(pcis->BtoNmap,IS_GTOLM_DROP,n_dir,idxs,&nmap,deluxe_ctx->idx_simple_B+n_com);CHKERRQ(ierr);
347       if (nmap != n_dir) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error when mapping simply scaled dofs (sub_schurs->is_dir)! %D != %D",nmap,n_dir);
348       ierr = ISRestoreIndices(sub_schurs->is_dir,&idxs);CHKERRQ(ierr);
349     }
350     ierr = PetscSortInt(deluxe_ctx->n_simple,deluxe_ctx->idx_simple_B);CHKERRQ(ierr);
351   } else {
352     deluxe_ctx->n_simple = 0;
353     deluxe_ctx->idx_simple_B = 0;
354   }
355   PetscFunctionReturn(0);
356 }
357 
358 #undef __FUNCT__
359 #define __FUNCT__ "PCBDDCScalingSetUp_Deluxe_Private"
360 static PetscErrorCode PCBDDCScalingSetUp_Deluxe_Private(PC pc)
361 {
362   PC_BDDC                *pcbddc=(PC_BDDC*)pc->data;
363   PCBDDCDeluxeScaling    deluxe_ctx=pcbddc->deluxe_ctx;
364   PCBDDCSubSchurs        sub_schurs = pcbddc->sub_schurs;
365   PetscErrorCode         ierr;
366 
367   PetscFunctionBegin;
368   if (!sub_schurs->n_subs) {
369     PetscFunctionReturn(0);
370   }
371 
372   /* Create work vectors for sequential part of deluxe */
373   ierr = MatCreateVecs(sub_schurs->S_Ej_all,&deluxe_ctx->seq_work1,&deluxe_ctx->seq_work2);CHKERRQ(ierr);
374 
375   /* Compute deluxe sequential scatter */
376   if (sub_schurs->reuse_mumps && !sub_schurs->is_dir) {
377     PCBDDCReuseMumps reuse_mumps = sub_schurs->reuse_mumps;
378     ierr = PetscObjectReference((PetscObject)reuse_mumps->correction_scatter_B);CHKERRQ(ierr);
379     deluxe_ctx->seq_scctx = reuse_mumps->correction_scatter_B;
380   } else {
381     ierr = VecScatterCreate(pcbddc->work_scaling,sub_schurs->is_Ej_all,deluxe_ctx->seq_work1,NULL,&deluxe_ctx->seq_scctx);CHKERRQ(ierr);
382   }
383 
384   /* Create Mat object for deluxe scaling */
385   ierr = PetscObjectReference((PetscObject)sub_schurs->S_Ej_all);CHKERRQ(ierr);
386   deluxe_ctx->seq_mat = sub_schurs->S_Ej_all;
387   if (sub_schurs->sum_S_Ej_all) { /* if this matrix is present, then we need to create the KSP object to invert it */
388     PC               pc_temp;
389     MatSolverPackage solver=NULL;
390     char             ksp_prefix[256];
391     size_t           len;
392 
393     ierr = KSPCreate(PETSC_COMM_SELF,&deluxe_ctx->seq_ksp);CHKERRQ(ierr);
394     ierr = KSPSetOperators(deluxe_ctx->seq_ksp,sub_schurs->sum_S_Ej_all,sub_schurs->sum_S_Ej_all);CHKERRQ(ierr);
395     ierr = KSPSetType(deluxe_ctx->seq_ksp,KSPPREONLY);CHKERRQ(ierr);
396     ierr = KSPGetPC(deluxe_ctx->seq_ksp,&pc_temp);CHKERRQ(ierr);
397     ierr = PCSetType(pc_temp,PCLU);CHKERRQ(ierr);
398     ierr = KSPGetPC(pcbddc->ksp_D,&pc_temp);CHKERRQ(ierr);
399     ierr = PCFactorGetMatSolverPackage(pc_temp,(const MatSolverPackage*)&solver);CHKERRQ(ierr);
400     if (solver) {
401       PC     new_pc;
402       PCType type;
403 
404       ierr = PCGetType(pc_temp,&type);CHKERRQ(ierr);
405       ierr = KSPGetPC(deluxe_ctx->seq_ksp,&new_pc);CHKERRQ(ierr);
406       ierr = PCSetType(new_pc,type);CHKERRQ(ierr);
407       ierr = PCFactorSetMatSolverPackage(new_pc,solver);CHKERRQ(ierr);
408     }
409     ierr = PetscStrlen(((PetscObject)(pcbddc->ksp_D))->prefix,&len);CHKERRQ(ierr);
410     len -= 10; /* remove "dirichlet_" */
411     ierr = PetscStrncpy(ksp_prefix,((PetscObject)(pcbddc->ksp_D))->prefix,len+1);CHKERRQ(ierr);
412     ierr = PetscStrcat(ksp_prefix,"deluxe_");CHKERRQ(ierr);
413     ierr = KSPSetOptionsPrefix(deluxe_ctx->seq_ksp,ksp_prefix);CHKERRQ(ierr);
414     ierr = KSPSetFromOptions(deluxe_ctx->seq_ksp);CHKERRQ(ierr);
415     ierr = KSPSetUp(deluxe_ctx->seq_ksp);CHKERRQ(ierr);
416   }
417   PetscFunctionReturn(0);
418 }
419