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