xref: /petsc/src/ksp/pc/impls/bddc/bddcscalingbasic.c (revision a07ea27a0d0aa5839ea22e119d3de928ef75918f) 
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   IS                  dirIS;
315   PetscErrorCode      ierr;
316 
317   PetscFunctionBegin;
318   /* (TODO: reuse) throw away the solvers */
319   ierr = PCBDDCScalingReset_Deluxe_Solvers(deluxe_ctx);CHKERRQ(ierr);
320 
321   /* Compute data structures to solve sequential problems */
322   ierr = PCBDDCScalingSetUp_Deluxe_Private(pc);CHKERRQ(ierr);
323 
324   /* diagonal scaling on interface dofs not contained in cc */
325   dirIS = NULL;
326   ierr = PCBDDCGraphGetDirichletDofs(pcbddc->mat_graph,&dirIS);CHKERRQ(ierr);
327   if (sub_schurs->is_Ej_com || dirIS) {
328     PetscInt n_com,n_dir;
329     n_com = 0;
330     if (sub_schurs->is_Ej_com) {
331       ierr = ISGetLocalSize(sub_schurs->is_Ej_com,&n_com);CHKERRQ(ierr);
332     }
333     n_dir = 0;
334     if (dirIS) {
335       ierr = ISGetLocalSize(dirIS,&n_dir);CHKERRQ(ierr);
336     }
337     deluxe_ctx->n_simple = n_dir + n_com;
338     ierr = PetscMalloc1(deluxe_ctx->n_simple,&deluxe_ctx->idx_simple_B);CHKERRQ(ierr);
339     if (sub_schurs->is_Ej_com) {
340       PetscInt       nmap;
341       const PetscInt *idxs;
342 
343       ierr = ISGetIndices(sub_schurs->is_Ej_com,&idxs);CHKERRQ(ierr);
344       ierr = ISGlobalToLocalMappingApply(pcis->BtoNmap,IS_GTOLM_DROP,n_com,idxs,&nmap,deluxe_ctx->idx_simple_B);CHKERRQ(ierr);
345       if (nmap != n_com) {
346         SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error when mapping simply scaled dofs (is_Ej_com)! %d != %d",nmap,n_com);
347       }
348       ierr = ISRestoreIndices(sub_schurs->is_Ej_com,&idxs);CHKERRQ(ierr);
349     }
350     if (dirIS) {
351       PetscInt       nmap;
352       const PetscInt *idxs;
353 
354       ierr = ISGetIndices(dirIS,&idxs);CHKERRQ(ierr);
355       ierr = ISGlobalToLocalMappingApply(pcis->BtoNmap,IS_GTOLM_DROP,n_dir,idxs,&nmap,deluxe_ctx->idx_simple_B+n_com);CHKERRQ(ierr);
356       ierr = ISRestoreIndices(dirIS,&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   ierr = ISDestroy(&dirIS);CHKERRQ(ierr);
364   PetscFunctionReturn(0);
365 }
366 
367 #undef __FUNCT__
368 #define __FUNCT__ "PCBDDCScalingSetUp_Deluxe_Private"
369 static PetscErrorCode PCBDDCScalingSetUp_Deluxe_Private(PC pc)
370 {
371   PC_BDDC                *pcbddc=(PC_BDDC*)pc->data;
372   PCBDDCDeluxeScaling    deluxe_ctx=pcbddc->deluxe_ctx;
373   PCBDDCSubSchurs        sub_schurs = pcbddc->sub_schurs;
374   PetscErrorCode         ierr;
375 
376   PetscFunctionBegin;
377   if (!sub_schurs->n_subs) {
378     PetscFunctionReturn(0);
379   }
380 
381   /* Create work vectors for sequential part of deluxe */
382   ierr = MatCreateVecs(sub_schurs->S_Ej_all,&deluxe_ctx->seq_work1,&deluxe_ctx->seq_work2);CHKERRQ(ierr);
383 
384   /* Compute deluxe sequential scatter */
385   ierr = VecScatterCreate(pcbddc->work_scaling,sub_schurs->is_Ej_all,deluxe_ctx->seq_work1,NULL,&deluxe_ctx->seq_scctx);CHKERRQ(ierr);
386 
387   /* Create Mat object for deluxe scaling */
388   ierr = PetscObjectReference((PetscObject)sub_schurs->S_Ej_all);CHKERRQ(ierr);
389   deluxe_ctx->seq_mat = sub_schurs->S_Ej_all;
390   if (sub_schurs->sum_S_Ej_all) { /* if this matrix is present, then we need to create the KSP object to invert it */
391     PC               pc_temp;
392     MatSolverPackage solver=NULL;
393     char             ksp_prefix[256];
394     size_t           len;
395 
396     ierr = KSPCreate(PETSC_COMM_SELF,&deluxe_ctx->seq_ksp);CHKERRQ(ierr);
397     ierr = KSPSetOperators(deluxe_ctx->seq_ksp,sub_schurs->sum_S_Ej_all,sub_schurs->sum_S_Ej_all);CHKERRQ(ierr);
398     ierr = KSPSetType(deluxe_ctx->seq_ksp,KSPPREONLY);CHKERRQ(ierr);
399     ierr = KSPGetPC(deluxe_ctx->seq_ksp,&pc_temp);CHKERRQ(ierr);
400     ierr = PCSetType(pc_temp,PCLU);CHKERRQ(ierr);
401     ierr = KSPGetPC(pcbddc->ksp_D,&pc_temp);CHKERRQ(ierr);
402     ierr = PCFactorGetMatSolverPackage(pc_temp,(const MatSolverPackage*)&solver);CHKERRQ(ierr);
403     if (solver) {
404       PC     new_pc;
405       PCType type;
406 
407       ierr = PCGetType(pc_temp,&type);CHKERRQ(ierr);
408       ierr = KSPGetPC(deluxe_ctx->seq_ksp,&new_pc);CHKERRQ(ierr);
409       ierr = PCSetType(new_pc,type);CHKERRQ(ierr);
410       ierr = PCFactorSetMatSolverPackage(new_pc,solver);CHKERRQ(ierr);
411     }
412     ierr = PetscStrlen(((PetscObject)(pcbddc->ksp_D))->prefix,&len);CHKERRQ(ierr);
413     len -= 10; /* remove "dirichlet_" */
414     ierr = PetscStrncpy(ksp_prefix,((PetscObject)(pcbddc->ksp_D))->prefix,len+1);CHKERRQ(ierr);
415     ierr = PetscStrcat(ksp_prefix,"deluxe_");CHKERRQ(ierr);
416     ierr = KSPSetOptionsPrefix(deluxe_ctx->seq_ksp,ksp_prefix);CHKERRQ(ierr);
417     ierr = KSPSetFromOptions(deluxe_ctx->seq_ksp);CHKERRQ(ierr);
418     ierr = KSPSetUp(deluxe_ctx->seq_ksp);CHKERRQ(ierr);
419   }
420   PetscFunctionReturn(0);
421 }
422