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