xref: /petsc/src/ksp/pc/impls/bddc/bddcscalingbasic.c (revision 310ea0fa8882f7442e116bd7f603f4ad71a130dc) 
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_Par(PC,PetscInt,PetscInt,PetscInt[],PetscInt[]);
9 static PetscErrorCode PCBDDCScalingSetUp_Deluxe_Seq(PC,PetscInt,PetscInt,PetscInt[],PetscInt[]);
10 static PetscErrorCode PCBDDCScalingReset_Deluxe_Solvers(PCBDDCDeluxeScaling);
11 
12 #undef __FUNCT__
13 #define __FUNCT__ "PCBDDCScalingExtension_Basic"
14 static PetscErrorCode PCBDDCScalingExtension_Basic(PC pc, Vec local_interface_vector, Vec global_vector)
15 {
16   PC_IS* pcis = (PC_IS*)pc->data;
17   PC_BDDC* pcbddc = (PC_BDDC*)pc->data;
18   PetscErrorCode ierr;
19 
20   PetscFunctionBegin;
21   /* Apply partition of unity */
22   ierr = VecPointwiseMult(pcbddc->work_scaling,pcis->D,local_interface_vector);CHKERRQ(ierr);
23   ierr = VecSet(global_vector,0.0);CHKERRQ(ierr);
24   ierr = VecScatterBegin(pcis->global_to_B,pcbddc->work_scaling,global_vector,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
25   ierr = VecScatterEnd(pcis->global_to_B,pcbddc->work_scaling,global_vector,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
26   PetscFunctionReturn(0);
27 }
28 
29 #undef __FUNCT__
30 #define __FUNCT__ "PCBDDCScalingExtension_Deluxe"
31 static PetscErrorCode PCBDDCScalingExtension_Deluxe(PC pc, Vec x, Vec y)
32 {
33   PC_IS*              pcis=(PC_IS*)pc->data;
34   PC_BDDC*            pcbddc=(PC_BDDC*)pc->data;
35   PCBDDCDeluxeScaling deluxe_ctx = pcbddc->deluxe_ctx;
36   PCBDDCSubSchurs     sub_schurs = deluxe_ctx->sub_schurs;
37   PetscInt            i;
38   PetscErrorCode      ierr;
39 
40   /* TODO CHECK STUFF RELATED WITH FAKE WORK */
41   PetscFunctionBegin;
42   ierr = VecSet(pcbddc->work_scaling,0.0);CHKERRQ(ierr); /* needed by the fake work below */
43   if (deluxe_ctx->n_simple) {
44     /* scale deluxe vertices using diagonal scaling */
45     PetscScalar *array_x,*array_D,*array;
46     ierr = VecGetArray(x,&array_x);CHKERRQ(ierr);
47     ierr = VecGetArray(pcis->D,&array_D);CHKERRQ(ierr);
48     ierr = VecGetArray(pcbddc->work_scaling,&array);CHKERRQ(ierr);
49     for (i=0;i<deluxe_ctx->n_simple;i++) {
50       array[deluxe_ctx->idx_simple_B[i]] = array_x[deluxe_ctx->idx_simple_B[i]]*array_D[deluxe_ctx->idx_simple_B[i]];
51     }
52     ierr = VecRestoreArray(pcbddc->work_scaling,&array);CHKERRQ(ierr);
53     ierr = VecRestoreArray(pcis->D,&array_D);CHKERRQ(ierr);
54     ierr = VecRestoreArray(x,&array_x);CHKERRQ(ierr);
55   }
56   /* sequential part : all problems and Schur applications collapsed into a single matrix vector multiplication and ksp solution */
57   if (deluxe_ctx->seq_mat) {
58     ierr = VecScatterBegin(deluxe_ctx->seq_scctx,x,deluxe_ctx->seq_work1,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
59     ierr = VecScatterEnd(deluxe_ctx->seq_scctx,x,deluxe_ctx->seq_work1,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
60     ierr = MatMult(deluxe_ctx->seq_mat,deluxe_ctx->seq_work1,deluxe_ctx->seq_work2);CHKERRQ(ierr);
61     ierr = KSPSolve(deluxe_ctx->seq_ksp,deluxe_ctx->seq_work2,deluxe_ctx->seq_work1);CHKERRQ(ierr);
62     /* fake work due to final ADD VALUES and vertices scaling needed? TODO: check it */
63     ierr = VecScatterBegin(deluxe_ctx->seq_scctx,deluxe_ctx->seq_work1,pcbddc->work_scaling,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
64     ierr = VecScatterEnd(deluxe_ctx->seq_scctx,deluxe_ctx->seq_work1,pcbddc->work_scaling,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
65   }
66   /* parallel part */
67   for (i=0;i<deluxe_ctx->par_colors;i++) {
68     if (deluxe_ctx->par_ksp[i]) {
69       PetscMPIInt color_rank;
70       PetscInt    subidx = deluxe_ctx->par_col2sub[i];
71       /* restrict on subset */
72       ierr = VecScatterBegin(deluxe_ctx->par_scctx_s[i],x,sub_schurs->work1[subidx],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
73       ierr = VecScatterEnd(deluxe_ctx->par_scctx_s[i],x,sub_schurs->work1[subidx],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
74       /* S_Ej */
75       ierr = MatMult(sub_schurs->S_Ej[subidx],sub_schurs->work1[subidx],sub_schurs->work2[subidx]);CHKERRQ(ierr);
76       /* (\sum_j S_Ej)^-1 */
77       ierr = VecSet(deluxe_ctx->par_vec[i],0.0);CHKERRQ(ierr);
78       ierr = VecScatterBegin(deluxe_ctx->par_scctx_p[i],sub_schurs->work2[subidx],deluxe_ctx->par_vec[i],ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
79       ierr = VecScatterEnd(deluxe_ctx->par_scctx_p[i],sub_schurs->work2[subidx],deluxe_ctx->par_vec[i],ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
80       ierr = KSPSolve(deluxe_ctx->par_ksp[i],deluxe_ctx->par_vec[i],deluxe_ctx->par_vec[i]);CHKERRQ(ierr);
81       ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)deluxe_ctx->par_ksp[i]),&color_rank);CHKERRQ(ierr);
82       /* get back solution on subset */
83       ierr = VecScatterBegin(deluxe_ctx->par_scctx_p[i],deluxe_ctx->par_vec[i],sub_schurs->work1[subidx],INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
84       ierr = VecScatterEnd(deluxe_ctx->par_scctx_p[i],deluxe_ctx->par_vec[i],sub_schurs->work1[subidx],INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
85       if (!color_rank) { /* only the master process in coloured comm copies the computed values */
86         ierr = VecScatterBegin(deluxe_ctx->par_scctx_s[i],sub_schurs->work1[subidx],pcbddc->work_scaling,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
87         ierr = VecScatterEnd(deluxe_ctx->par_scctx_s[i],sub_schurs->work1[subidx],pcbddc->work_scaling,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
88       }
89     }
90   }
91   /* put local boundary part in global vector */
92   ierr = VecSet(y,0.0);CHKERRQ(ierr);
93   ierr = VecScatterBegin(pcis->global_to_B,pcbddc->work_scaling,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
94   ierr = VecScatterEnd(pcis->global_to_B,pcbddc->work_scaling,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
95   PetscFunctionReturn(0);
96 }
97 
98 #undef __FUNCT__
99 #define __FUNCT__ "PCBDDCScalingExtension"
100 PetscErrorCode PCBDDCScalingExtension(PC pc, Vec local_interface_vector, Vec global_vector)
101 {
102   PC_BDDC *pcbddc=(PC_BDDC*)pc->data;
103   PetscErrorCode ierr;
104 
105   PetscFunctionBegin;
106   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
107   PetscValidHeaderSpecific(local_interface_vector,VEC_CLASSID,2);
108   PetscValidHeaderSpecific(global_vector,VEC_CLASSID,3);
109   if (local_interface_vector == pcbddc->work_scaling) {
110     SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Local vector cannot be pcbddc->work_scaling!\n");
111   }
112   ierr = PetscTryMethod(pc,"PCBDDCScalingExtension_C",(PC,Vec,Vec),(pc,local_interface_vector,global_vector));CHKERRQ(ierr);
113   PetscFunctionReturn(0);
114 }
115 
116 #undef __FUNCT__
117 #define __FUNCT__ "PCBDDCScalingRestriction_Basic"
118 static PetscErrorCode PCBDDCScalingRestriction_Basic(PC pc, Vec global_vector, Vec local_interface_vector)
119 {
120   PetscErrorCode ierr;
121   PC_IS* pcis = (PC_IS*)pc->data;
122 
123   PetscFunctionBegin;
124   ierr = VecScatterBegin(pcis->global_to_B,global_vector,local_interface_vector,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
125   ierr = VecScatterEnd(pcis->global_to_B,global_vector,local_interface_vector,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
126   /* Apply partition of unity */
127   ierr = VecPointwiseMult(local_interface_vector,pcis->D,local_interface_vector);CHKERRQ(ierr);
128   PetscFunctionReturn(0);
129 }
130 
131 #undef __FUNCT__
132 #define __FUNCT__ "PCBDDCScalingRestriction_Deluxe"
133 static PetscErrorCode PCBDDCScalingRestriction_Deluxe(PC pc, Vec x, Vec y)
134 {
135   PC_IS*              pcis=(PC_IS*)pc->data;
136   PC_BDDC*            pcbddc=(PC_BDDC*)pc->data;
137   PCBDDCDeluxeScaling deluxe_ctx = pcbddc->deluxe_ctx;
138   PCBDDCSubSchurs     sub_schurs = deluxe_ctx->sub_schurs;
139   PetscInt            i;
140   PetscErrorCode      ierr;
141 
142   PetscFunctionBegin;
143   /* get local boundary part of global vector */
144   ierr = VecScatterBegin(pcis->global_to_B,x,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
145   ierr = VecScatterEnd(pcis->global_to_B,x,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
146   if (deluxe_ctx->n_simple) {
147     /* scale deluxe vertices using diagonal scaling */
148     PetscScalar *array_y,*array_D;
149     ierr = VecGetArray(y,&array_y);CHKERRQ(ierr);
150     ierr = VecGetArray(pcis->D,&array_D);CHKERRQ(ierr);
151     for (i=0;i<deluxe_ctx->n_simple;i++) {
152       array_y[deluxe_ctx->idx_simple_B[i]] *= array_D[deluxe_ctx->idx_simple_B[i]];
153     }
154     ierr = VecRestoreArray(pcis->D,&array_D);CHKERRQ(ierr);
155     ierr = VecRestoreArray(y,&array_y);CHKERRQ(ierr);
156   }
157   /* sequential part : all problems and Schur applications collapsed into a single matrix vector multiplication and ksp solution */
158   if (deluxe_ctx->seq_mat) {
159     ierr = VecScatterBegin(deluxe_ctx->seq_scctx,y,deluxe_ctx->seq_work1,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
160     ierr = VecScatterEnd(deluxe_ctx->seq_scctx,y,deluxe_ctx->seq_work1,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
161     ierr = KSPSolveTranspose(deluxe_ctx->seq_ksp,deluxe_ctx->seq_work1,deluxe_ctx->seq_work2);CHKERRQ(ierr);
162     ierr = MatMultTranspose(deluxe_ctx->seq_mat,deluxe_ctx->seq_work2,deluxe_ctx->seq_work1);CHKERRQ(ierr);
163     ierr = VecScatterBegin(deluxe_ctx->seq_scctx,deluxe_ctx->seq_work1,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
164     ierr = VecScatterEnd(deluxe_ctx->seq_scctx,deluxe_ctx->seq_work1,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
165   }
166   /* parallel part */
167   for (i=0;i<deluxe_ctx->par_colors;i++) {
168     if (deluxe_ctx->par_ksp[i]) {
169       PetscInt subidx = deluxe_ctx->par_col2sub[i];
170       /* restrict on subset */
171       ierr = VecScatterBegin(deluxe_ctx->par_scctx_s[i],y,sub_schurs->work1[subidx],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
172       ierr = VecScatterEnd(deluxe_ctx->par_scctx_s[i],y,sub_schurs->work1[subidx],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
173       /* (\sum_j S_Ej)^-T */
174       ierr = VecSet(deluxe_ctx->par_vec[i],0.0);CHKERRQ(ierr);
175       ierr = VecScatterBegin(deluxe_ctx->par_scctx_p[i],sub_schurs->work1[subidx],deluxe_ctx->par_vec[i],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
176       ierr = VecScatterEnd(deluxe_ctx->par_scctx_p[i],sub_schurs->work1[subidx],deluxe_ctx->par_vec[i],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
177       ierr = KSPSolveTranspose(deluxe_ctx->par_ksp[i],deluxe_ctx->par_vec[i],deluxe_ctx->par_vec[i]);CHKERRQ(ierr);
178       ierr = VecScatterBegin(deluxe_ctx->par_scctx_p[i],deluxe_ctx->par_vec[i],sub_schurs->work1[subidx],INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
179       ierr = VecScatterEnd(deluxe_ctx->par_scctx_p[i],deluxe_ctx->par_vec[i],sub_schurs->work1[subidx],INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
180       /* S_Ej^T */
181       ierr = MatMultTranspose(sub_schurs->S_Ej[subidx],sub_schurs->work1[subidx],sub_schurs->work2[subidx]);CHKERRQ(ierr);
182       /* extend to boundary */
183       ierr = VecScatterBegin(deluxe_ctx->par_scctx_s[i],sub_schurs->work2[subidx],y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
184       ierr = VecScatterEnd(deluxe_ctx->par_scctx_s[i],sub_schurs->work2[subidx],y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
185     }
186   }
187   PetscFunctionReturn(0);
188 }
189 
190 #undef __FUNCT__
191 #define __FUNCT__ "PCBDDCScalingRestriction"
192 PetscErrorCode PCBDDCScalingRestriction(PC pc, Vec global_vector, Vec local_interface_vector)
193 {
194   PC_BDDC        *pcbddc=(PC_BDDC*)pc->data;
195   PetscErrorCode ierr;
196 
197   PetscFunctionBegin;
198   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
199   PetscValidHeaderSpecific(global_vector,VEC_CLASSID,2);
200   PetscValidHeaderSpecific(local_interface_vector,VEC_CLASSID,3);
201   if (local_interface_vector == pcbddc->work_scaling) {
202     SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Local vector should cannot be pcbddc->work_scaling!\n");
203   }
204   ierr = PetscTryMethod(pc,"PCBDDCScalingRestriction_C",(PC,Vec,Vec),(pc,global_vector,local_interface_vector));CHKERRQ(ierr);
205   PetscFunctionReturn(0);
206 }
207 
208 #undef __FUNCT__
209 #define __FUNCT__ "PCBDDCScalingSetUp"
210 PetscErrorCode PCBDDCScalingSetUp(PC pc)
211 {
212   PC_IS* pcis=(PC_IS*)pc->data;
213   PC_BDDC* pcbddc=(PC_BDDC*)pc->data;
214   PetscErrorCode ierr;
215 
216   PetscFunctionBegin;
217   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
218   /* create work vector for the operator */
219   ierr = VecDestroy(&pcbddc->work_scaling);CHKERRQ(ierr);
220   ierr = VecDuplicate(pcis->vec1_B,&pcbddc->work_scaling);CHKERRQ(ierr);
221   /* always rebuild pcis->D */
222   if (pcis->use_stiffness_scaling) {
223     ierr = MatGetDiagonal(pcbddc->local_mat,pcis->vec1_N);CHKERRQ(ierr);
224     ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
225     ierr = VecScatterEnd(pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
226   }
227   ierr = VecCopy(pcis->D,pcis->vec1_B);CHKERRQ(ierr);
228   ierr = VecSet(pcis->vec1_global,0.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 = VecScatterBegin(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
232   ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
233   ierr = VecPointwiseDivide(pcis->D,pcis->D,pcis->vec1_B);CHKERRQ(ierr);
234   /* now setup */
235   if (pcbddc->use_deluxe_scaling) {
236     if (!pcbddc->deluxe_ctx) {
237       ierr = PCBDDCScalingCreate_Deluxe(pc);CHKERRQ(ierr);
238     }
239     ierr = PCBDDCScalingSetUp_Deluxe(pc);CHKERRQ(ierr);
240     ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingRestriction_C",PCBDDCScalingRestriction_Deluxe);CHKERRQ(ierr);
241     ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingExtension_C",PCBDDCScalingExtension_Deluxe);CHKERRQ(ierr);
242   } else {
243     ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingRestriction_C",PCBDDCScalingRestriction_Basic);CHKERRQ(ierr);
244     ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingExtension_C",PCBDDCScalingExtension_Basic);CHKERRQ(ierr);
245   }
246 
247   /* test */
248   if (pcbddc->dbg_flag) {
249     Vec         vec2_global;
250     PetscViewer viewer=pcbddc->dbg_viewer;
251     PetscReal   error;
252 
253     /* extension -> from local to parallel */
254     ierr = VecSet(pcis->vec1_global,0.0);CHKERRQ(ierr);
255     ierr = VecSetRandom(pcis->vec1_B,NULL);CHKERRQ(ierr);
256     ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
257     ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
258     ierr = VecDuplicate(pcis->vec1_global,&vec2_global);CHKERRQ(ierr);
259     ierr = VecCopy(pcis->vec1_global,vec2_global);CHKERRQ(ierr);
260 
261     ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
262     ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
263     ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,pcis->vec1_global);CHKERRQ(ierr);
264     ierr = VecAXPY(pcis->vec1_global,-1.0,vec2_global);CHKERRQ(ierr);
265     ierr = VecNorm(pcis->vec1_global,NORM_INFINITY,&error);CHKERRQ(ierr);
266     ierr = PetscViewerASCIIPrintf(viewer,"Error scaling extension %1.14e\n",error);CHKERRQ(ierr);
267     if (error>1.e-8 && pcbddc->dbg_flag>1) {
268       ierr = VecView(pcis->vec1_global,viewer);CHKERRQ(ierr);
269     }
270     ierr = VecDestroy(&vec2_global);CHKERRQ(ierr);
271 
272     /* restriction -> from parallel to local */
273     ierr = VecSet(pcis->vec1_global,0.0);CHKERRQ(ierr);
274     ierr = VecSetRandom(pcis->vec1_B,NULL);CHKERRQ(ierr);
275     ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
276     ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
277 
278     ierr = PCBDDCScalingRestriction(pc,pcis->vec1_global,pcis->vec1_B);CHKERRQ(ierr);
279     ierr = VecScale(pcis->vec1_B,-1.0);CHKERRQ(ierr);
280     ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
281     ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
282     ierr = VecNorm(pcis->vec1_global,NORM_INFINITY,&error);CHKERRQ(ierr);
283     ierr = PetscViewerASCIIPrintf(viewer,"Error scaling restriction %1.14e\n",error);CHKERRQ(ierr);
284     if (error>1.e-8 && pcbddc->dbg_flag>1) {
285       ierr = VecView(pcis->vec1_global,viewer);CHKERRQ(ierr);
286     }
287   }
288   PetscFunctionReturn(0);
289 }
290 
291 #undef __FUNCT__
292 #define __FUNCT__ "PCBDDCScalingDestroy"
293 PetscErrorCode PCBDDCScalingDestroy(PC pc)
294 {
295   PC_BDDC* pcbddc=(PC_BDDC*)pc->data;
296   PetscErrorCode ierr;
297 
298   PetscFunctionBegin;
299   if (pcbddc->deluxe_ctx) {
300     ierr = PCBDDCScalingDestroy_Deluxe(pc);CHKERRQ(ierr);
301   }
302   ierr = VecDestroy(&pcbddc->work_scaling);CHKERRQ(ierr);
303   /* remove functions */
304   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingRestriction_C",NULL);CHKERRQ(ierr);
305   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingExtension_C",NULL);CHKERRQ(ierr);
306   PetscFunctionReturn(0);
307 }
308 
309 #undef __FUNCT__
310 #define __FUNCT__ "PCBDDCScalingCreate_Deluxe"
311 static PetscErrorCode PCBDDCScalingCreate_Deluxe(PC pc)
312 {
313   PC_BDDC*            pcbddc=(PC_BDDC*)pc->data;
314   PCBDDCDeluxeScaling deluxe_ctx;
315   PetscErrorCode      ierr;
316 
317   PetscFunctionBegin;
318   ierr = PetscNew(&deluxe_ctx);CHKERRQ(ierr);
319   ierr = PCBDDCSubSchursCreate(&deluxe_ctx->sub_schurs);CHKERRQ(ierr);
320   pcbddc->deluxe_ctx = deluxe_ctx;
321   PetscFunctionReturn(0);
322 }
323 
324 #undef __FUNCT__
325 #define __FUNCT__ "PCBDDCScalingDestroy_Deluxe"
326 static PetscErrorCode PCBDDCScalingDestroy_Deluxe(PC pc)
327 {
328   PC_BDDC*            pcbddc=(PC_BDDC*)pc->data;
329   PetscErrorCode      ierr;
330 
331   PetscFunctionBegin;
332   ierr = PCBDDCScalingReset_Deluxe_Solvers(pcbddc->deluxe_ctx);CHKERRQ(ierr);
333   ierr = PCBDDCSubSchursDestroy(&(pcbddc->deluxe_ctx->sub_schurs));CHKERRQ(ierr);
334   ierr = PetscFree(pcbddc->deluxe_ctx);CHKERRQ(ierr);
335   PetscFunctionReturn(0);
336 }
337 
338 #undef __FUNCT__
339 #define __FUNCT__ "PCBDDCScalingReset_Deluxe_Solvers"
340 static PetscErrorCode PCBDDCScalingReset_Deluxe_Solvers(PCBDDCDeluxeScaling deluxe_ctx)
341 {
342   PetscErrorCode      ierr;
343 
344   PetscFunctionBegin;
345   ierr = PetscFree(deluxe_ctx->idx_simple_B);CHKERRQ(ierr);
346   deluxe_ctx->n_simple = 0;
347   if (deluxe_ctx->seq_mat) {
348     ierr = VecScatterDestroy(&deluxe_ctx->seq_scctx);CHKERRQ(ierr);
349     ierr = VecDestroy(&deluxe_ctx->seq_work1);CHKERRQ(ierr);
350     ierr = VecDestroy(&deluxe_ctx->seq_work2);CHKERRQ(ierr);
351     ierr = MatDestroy(&deluxe_ctx->seq_mat);CHKERRQ(ierr);
352     ierr = KSPDestroy(&deluxe_ctx->seq_ksp);CHKERRQ(ierr);
353   }
354   if (deluxe_ctx->par_colors) {
355     PetscInt i;
356     for (i=0;i<deluxe_ctx->par_colors;i++) {
357       ierr = VecScatterDestroy(&deluxe_ctx->par_scctx_s[i]);CHKERRQ(ierr);
358       ierr = VecScatterDestroy(&deluxe_ctx->par_scctx_p[i]);CHKERRQ(ierr);
359       ierr = VecDestroy(&deluxe_ctx->par_vec[i]);CHKERRQ(ierr);
360       ierr = KSPDestroy(&deluxe_ctx->par_ksp[i]);CHKERRQ(ierr);
361     }
362     ierr = PetscFree5(deluxe_ctx->par_ksp,
363                       deluxe_ctx->par_scctx_s,
364                       deluxe_ctx->par_scctx_p,
365                       deluxe_ctx->par_vec,
366                       deluxe_ctx->par_col2sub);CHKERRQ(ierr);
367   }
368   deluxe_ctx->par_colors = 0;
369   PetscFunctionReturn(0);
370 }
371 
372 #undef __FUNCT__
373 #define __FUNCT__ "PCBDDCScalingSetUp_Deluxe"
374 static PetscErrorCode PCBDDCScalingSetUp_Deluxe(PC pc)
375 {
376   PC_IS               *pcis=(PC_IS*)pc->data;
377   PC_BDDC             *pcbddc=(PC_BDDC*)pc->data;
378   PCBDDCDeluxeScaling deluxe_ctx=pcbddc->deluxe_ctx;
379   PCBDDCSubSchurs     sub_schurs=deluxe_ctx->sub_schurs;
380   PCBDDCGraph         graph;
381   IS                  *faces,*edges,*all_cc;
382   PetscBT             bitmask;
383   PetscInt            *index_sequential,*index_parallel;
384   PetscInt            *auxlocal_sequential,*auxlocal_parallel;
385   PetscInt            *auxglobal_sequential,*auxglobal_parallel;
386   PetscInt            *auxmapping,*idxs;
387   PetscInt            i,max_subset_size;
388   PetscInt            n_sequential_problems,n_local_sequential_problems,n_parallel_problems,n_local_parallel_problems;
389   PetscInt            n_faces,n_edges,n_all_cc;
390   PetscErrorCode      ierr;
391 
392   PetscFunctionBegin;
393   /* throw away the solvers */
394   ierr = PCBDDCScalingReset_Deluxe_Solvers(deluxe_ctx);CHKERRQ(ierr);
395 
396   /* attach interface graph for determining subsets */
397   if (pcbddc->deluxe_rebuild) { /* in case rebuild has been requested, it uses a graph generated only by the neighbouring information */
398     PetscInt *idx_V_N;
399     IS       verticesIS;
400     ierr = PCBDDCGetPrimalVerticesLocalIdx(pc,&i,&idx_V_N);CHKERRQ(ierr);
401     ierr = ISCreateGeneral(PETSC_COMM_SELF,i,idx_V_N,PETSC_OWN_POINTER,&verticesIS);CHKERRQ(ierr);
402     ierr = PCBDDCGraphCreate(&graph);CHKERRQ(ierr);
403     ierr = PCBDDCGraphInit(graph,pcbddc->mat_graph->l2gmap);CHKERRQ(ierr);
404     ierr = PCBDDCGraphSetUp(graph,0,NULL,pcbddc->DirichletBoundariesLocal,0,NULL,verticesIS);CHKERRQ(ierr);
405     ierr = PCBDDCGraphComputeConnectedComponents(graph);CHKERRQ(ierr);
406     ierr = ISDestroy(&verticesIS);CHKERRQ(ierr);
407 /*
408     if (pcbddc->dbg_flag) {
409       ierr = PCBDDCGraphASCIIView(graph,pcbddc->dbg_flag,pcbddc->dbg_viewer);CHKERRQ(ierr);
410     }
411 */
412   } else {
413     graph = pcbddc->mat_graph;
414   }
415 
416   /* get index sets for faces and edges */
417   ierr = PCBDDCGraphGetCandidatesIS(graph,&n_faces,&faces,&n_edges,&edges,NULL);CHKERRQ(ierr);
418   n_all_cc = n_faces+n_edges;
419   ierr = PetscMalloc1(n_all_cc,&all_cc);CHKERRQ(ierr);
420   for (i=0;i<n_faces;i++) {
421     all_cc[i] = faces[i];
422   }
423   for (i=0;i<n_edges;i++) {
424     all_cc[n_faces+i] = edges[i];
425   }
426   ierr = PetscFree(faces);CHKERRQ(ierr);
427   ierr = PetscFree(edges);CHKERRQ(ierr);
428 
429   /* map interface's subsets */
430   max_subset_size = 0;
431   for (i=0;i<n_all_cc;i++) {
432     PetscInt subset_size;
433     ierr = ISGetLocalSize(all_cc[i],&subset_size);CHKERRQ(ierr);
434     max_subset_size = PetscMax(max_subset_size,subset_size);
435   }
436   ierr = PetscMalloc5(max_subset_size,&auxmapping,
437                       graph->ncc,&auxlocal_sequential,
438                       graph->ncc,&auxlocal_parallel,
439                       graph->ncc,&index_sequential,
440                       graph->ncc,&index_parallel);CHKERRQ(ierr);
441 
442   /* if threshold is negative, uses all sequential problems */
443   if (pcbddc->deluxe_threshold < 0) pcbddc->deluxe_threshold = max_subset_size;
444 
445   /* workspace */
446   ierr = PetscBTCreate(pcis->n,&bitmask);CHKERRQ(ierr);
447   ierr = ISGetIndices(pcis->is_I_local,(const PetscInt**)&idxs);CHKERRQ(ierr);
448   for (i=0;i<pcis->n-pcis->n_B;i++) {
449     ierr = PetscBTSet(bitmask,idxs[i]);CHKERRQ(ierr);
450   }
451   ierr = ISRestoreIndices(pcis->is_I_local,(const PetscInt**)&idxs);CHKERRQ(ierr);
452 
453   /* determine which problem has to be solved in parallel or sequentially */
454   n_local_sequential_problems = 0;
455   n_local_parallel_problems = 0;
456   for (i=0;i<n_all_cc;i++) {
457     PetscInt subset_size,j,min_loc = 0;
458 
459     ierr = ISGetLocalSize(all_cc[i],&subset_size);CHKERRQ(ierr);
460     ierr = ISGetIndices(all_cc[i],(const PetscInt**)&idxs);CHKERRQ(ierr);
461     for (j=0;j<subset_size;j++) {
462       ierr = PetscBTSet(bitmask,idxs[j]);CHKERRQ(ierr);
463     }
464     ierr = ISLocalToGlobalMappingApply(graph->l2gmap,subset_size,idxs,auxmapping);CHKERRQ(ierr);
465     for (j=1;j<subset_size;j++) {
466       if (auxmapping[j]<auxmapping[min_loc]) {
467         min_loc = j;
468       }
469     }
470     if (subset_size > pcbddc->deluxe_threshold) {
471       index_parallel[n_local_parallel_problems] = i;
472       auxlocal_parallel[n_local_parallel_problems] = idxs[min_loc];
473       n_local_parallel_problems++;
474     } else {
475       index_sequential[n_local_sequential_problems] = i;
476       auxlocal_sequential[n_local_sequential_problems] = idxs[min_loc];
477       n_local_sequential_problems++;
478     }
479     ierr = ISRestoreIndices(all_cc[i],(const PetscInt**)&idxs);CHKERRQ(ierr);
480   }
481 
482   /* diagonal scaling on interface dofs not contained in cc */
483   deluxe_ctx->n_simple = 0;
484   for (i=0;i<pcis->n;i++) {
485     if (!PetscBTLookup(bitmask,i)) {
486       deluxe_ctx->n_simple++;
487     }
488   }
489   ierr = PetscMalloc1(deluxe_ctx->n_simple,&deluxe_ctx->idx_simple_B);CHKERRQ(ierr);
490   deluxe_ctx->n_simple = 0;
491   for (i=0;i<pcis->n;i++) {
492     if (!PetscBTLookup(bitmask,i)) {
493       deluxe_ctx->idx_simple_B[deluxe_ctx->n_simple++] = i;
494     }
495   }
496   ierr = ISGlobalToLocalMappingApply(pcbddc->BtoNmap,IS_GTOLM_DROP,deluxe_ctx->n_simple,deluxe_ctx->idx_simple_B,&i,deluxe_ctx->idx_simple_B);CHKERRQ(ierr);
497   if (i != deluxe_ctx->n_simple) {
498     SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error when mapping simple scaling dofs! %d != %d",i,deluxe_ctx->n_simple);
499   }
500   ierr = PetscBTDestroy(&bitmask);CHKERRQ(ierr);
501 
502   /* SetUp local schur complements on subsets TODO better reuse procedure */
503   if (!sub_schurs->n_subs) {
504     Mat       S_j;
505     PetscBool free_used_adj;
506     PetscInt  *used_xadj,*used_adjncy;
507 
508     /* decide the adjacency to be used for determining internal problems for local schur on subsets */
509     free_used_adj = PETSC_FALSE;
510     if (pcbddc->deluxe_layers == -1) {
511       used_xadj = NULL;
512       used_adjncy = NULL;
513     } else {
514       if ((pcbddc->deluxe_use_useradj && pcbddc->mat_graph->xadj) || !pcbddc->deluxe_compute_rowadj) {
515         used_xadj = pcbddc->mat_graph->xadj;
516         used_adjncy = pcbddc->mat_graph->adjncy;
517       } else {
518         Mat            mat_adj;
519         PetscBool      flg_row=PETSC_TRUE;
520         const PetscInt *xadj,*adjncy;
521         PetscInt       nvtxs;
522 
523         ierr = MatConvert(pcbddc->local_mat,MATMPIADJ,MAT_INITIAL_MATRIX,&mat_adj);CHKERRQ(ierr);
524         ierr = MatGetRowIJ(mat_adj,0,PETSC_TRUE,PETSC_FALSE,&nvtxs,&xadj,&adjncy,&flg_row);CHKERRQ(ierr);
525         if (!flg_row) {
526           SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in MatGetRowIJ called in %s\n",__FUNCT__);
527         }
528         ierr = PetscMalloc2(nvtxs+1,&used_xadj,xadj[nvtxs],&used_adjncy);CHKERRQ(ierr);
529         ierr = PetscMemcpy(used_xadj,xadj,(nvtxs+1)*sizeof(*xadj));CHKERRQ(ierr);
530         ierr = PetscMemcpy(used_adjncy,adjncy,(xadj[nvtxs])*sizeof(*adjncy));CHKERRQ(ierr);
531         ierr = MatRestoreRowIJ(mat_adj,0,PETSC_TRUE,PETSC_FALSE,&nvtxs,&xadj,&adjncy,&flg_row);CHKERRQ(ierr);
532         if (!flg_row) {
533           SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in MatRestoreRowIJ called in %s\n",__FUNCT__);
534         }
535         ierr = MatDestroy(&mat_adj);CHKERRQ(ierr);
536         free_used_adj = PETSC_TRUE;
537       }
538     }
539 
540     /* Create Schur complement matrix */
541     ierr = MatCreateSchurComplement(pcis->A_II,pcis->A_II,pcis->A_IB,pcis->A_BI,pcis->A_BB,&S_j);CHKERRQ(ierr);
542     ierr = MatSchurComplementSetKSP(S_j,pcbddc->ksp_D);CHKERRQ(ierr);
543 
544     /* setup Schur complements on subsets */
545     ierr = PCBDDCSubSchursSetUp(sub_schurs,S_j,pcis->is_I_local,pcis->is_B_local,n_all_cc,all_cc,used_xadj,used_adjncy,pcbddc->deluxe_layers);CHKERRQ(ierr);
546     ierr = MatDestroy(&S_j);CHKERRQ(ierr);
547     /* free adjacency */
548     if (free_used_adj) {
549       ierr = PetscFree2(used_xadj,used_adjncy);CHKERRQ(ierr);
550     }
551   }
552   for (i=0;i<n_all_cc;i++) {
553     ierr = ISDestroy(&all_cc[i]);CHKERRQ(ierr);
554   }
555   ierr = PetscFree(all_cc);CHKERRQ(ierr);
556 
557   /* Number parallel problems */
558   auxglobal_parallel = 0;
559   ierr = PCBDDCSubsetNumbering(PetscObjectComm((PetscObject)pc),graph->l2gmap,n_local_parallel_problems,auxlocal_parallel,PETSC_NULL,&n_parallel_problems,&auxglobal_parallel);CHKERRQ(ierr);
560   if (pcbddc->dbg_flag) {
561     ierr = PetscViewerASCIIPrintf(pcbddc->dbg_viewer,"Deluxe global number of parallel subproblems: %d\n",n_parallel_problems);
562   }
563 
564   /* Compute data structures to solve parallel problems */
565   ierr = PCBDDCScalingSetUp_Deluxe_Par(pc,n_local_parallel_problems,n_parallel_problems,auxglobal_parallel,index_parallel);CHKERRQ(ierr);
566   ierr = PetscFree(auxglobal_parallel);CHKERRQ(ierr);
567 
568 
569   /* Number sequential problems */
570   auxglobal_sequential = 0;
571   ierr = PCBDDCSubsetNumbering(PetscObjectComm((PetscObject)pc),graph->l2gmap,n_local_sequential_problems,auxlocal_sequential,PETSC_NULL,&n_sequential_problems,&auxglobal_sequential);CHKERRQ(ierr);
572   if (pcbddc->dbg_flag) {
573     ierr = PetscViewerASCIIPrintf(pcbddc->dbg_viewer,"Deluxe global number of sequential subproblems: %d\n",n_sequential_problems);
574   }
575 
576   /* Compute data structures to solve sequential problems */
577   ierr = PCBDDCScalingSetUp_Deluxe_Seq(pc,n_local_sequential_problems,n_sequential_problems,auxglobal_sequential,index_sequential);CHKERRQ(ierr);
578   ierr = PetscFree(auxglobal_sequential);CHKERRQ(ierr);
579 
580   /* free workspace */
581   ierr = PetscFree5(auxmapping,auxlocal_sequential,auxlocal_parallel,index_sequential,index_parallel);CHKERRQ(ierr);
582 
583   /* free graph struct */
584   if (pcbddc->deluxe_rebuild) {
585     ierr = PCBDDCGraphDestroy(&graph);CHKERRQ(ierr);
586   }
587   PetscFunctionReturn(0);
588 }
589 
590 #undef __FUNCT__
591 #define __FUNCT__ "PCBDDCScalingSetUp_Deluxe_Par"
592 static PetscErrorCode PCBDDCScalingSetUp_Deluxe_Par(PC pc, PetscInt n_local_parallel_problems,PetscInt n_parallel_problems,PetscInt global_parallel[],PetscInt index_parallel[])
593 {
594   PC_BDDC                *pcbddc=(PC_BDDC*)pc->data;
595   PCBDDCDeluxeScaling    deluxe_ctx=pcbddc->deluxe_ctx;
596   /* coloring */
597   Mat                    parallel_problems;
598   MatColoring            coloring_obj;
599   ISColoring             coloring_parallel_problems;
600   IS                     *par_is_colors,*is_colors;
601   /* working stuff */
602   PetscInt               i,j;
603   PetscErrorCode         ierr;
604 
605   PetscFunctionBegin;
606   if (!n_parallel_problems) {
607     PetscFunctionReturn(0);
608   }
609   /* Color parallel subproblems */
610   ierr = MatCreate(PetscObjectComm((PetscObject)pc),&parallel_problems);CHKERRQ(ierr);
611   ierr = MatSetSizes(parallel_problems,PETSC_DECIDE,PETSC_DECIDE,n_parallel_problems,n_parallel_problems);CHKERRQ(ierr);
612   ierr = MatSetType(parallel_problems,MATAIJ);CHKERRQ(ierr);
613   ierr = MatSetUp(parallel_problems);CHKERRQ(ierr);
614   ierr = MatSetOption(parallel_problems,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr);
615   ierr = MatSetOption(parallel_problems,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr);
616   for (i=0;i<n_local_parallel_problems;i++) {
617     PetscInt row = global_parallel[i];
618     for (j=0;j<n_local_parallel_problems;j++) {
619       PetscInt col = global_parallel[j];
620       if (row != col) {
621         ierr = MatSetValue(parallel_problems,row,col,1.0,INSERT_VALUES);CHKERRQ(ierr);
622       }
623     }
624   }
625   ierr = MatAssemblyBegin(parallel_problems,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
626   ierr = MatAssemblyEnd(parallel_problems,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
627   if (pcbddc->dbg_flag > 1) {
628     ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr);
629     ierr = PetscViewerASCIIPrintf(pcbddc->dbg_viewer,"Adj matrix for deluxe parallel problems\n");CHKERRQ(ierr);
630     ierr = MatView(parallel_problems,pcbddc->dbg_viewer);CHKERRQ(ierr);
631   }
632   ierr = MatColoringCreate(parallel_problems,&coloring_obj);CHKERRQ(ierr);
633   ierr = MatColoringSetDistance(coloring_obj,1);CHKERRQ(ierr);
634   ierr = MatColoringSetType(coloring_obj,MATCOLORINGJP);CHKERRQ(ierr);
635   ierr = MatColoringApply(coloring_obj,&coloring_parallel_problems);CHKERRQ(ierr);
636   ierr = ISColoringGetIS(coloring_parallel_problems,&deluxe_ctx->par_colors,&par_is_colors);CHKERRQ(ierr);
637   if (pcbddc->dbg_flag) {
638     ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr);
639     ierr = PetscViewerASCIIPrintf(pcbddc->dbg_viewer,"Number of colors %d for parallel part of deluxe\n",deluxe_ctx->par_colors);CHKERRQ(ierr);
640   }
641 
642   /* all procs should know the color distribution */
643   ierr = PetscMalloc1(deluxe_ctx->par_colors,&is_colors);CHKERRQ(ierr);
644   for (i=0;i<deluxe_ctx->par_colors;i++) {
645     if (pcbddc->dbg_flag) {
646       ierr = PetscViewerASCIIPrintf(pcbddc->dbg_viewer,"Global problem indexes for color %d\n",i);CHKERRQ(ierr);
647       ierr = ISView(par_is_colors[i],pcbddc->dbg_viewer);CHKERRQ(ierr);
648       ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr);
649     }
650     ierr = ISAllGather(par_is_colors[i],&is_colors[i]);CHKERRQ(ierr);
651   }
652 
653   /* free unneeded objects */
654   ierr = ISColoringRestoreIS(coloring_parallel_problems,&par_is_colors);CHKERRQ(ierr);
655   ierr = ISColoringDestroy(&coloring_parallel_problems);CHKERRQ(ierr);
656   ierr = MatColoringDestroy(&coloring_obj);CHKERRQ(ierr);
657   ierr = MatDestroy(&parallel_problems);CHKERRQ(ierr);
658 
659   /* allocate deluxe arrays for parallel problems */
660   ierr = PetscMalloc5(deluxe_ctx->par_colors,&deluxe_ctx->par_ksp,
661                       deluxe_ctx->par_colors,&deluxe_ctx->par_scctx_s,
662                       deluxe_ctx->par_colors,&deluxe_ctx->par_scctx_p,
663                       deluxe_ctx->par_colors,&deluxe_ctx->par_vec,
664                       deluxe_ctx->par_colors,&deluxe_ctx->par_col2sub);CHKERRQ(ierr);
665 
666   /* cycle on colors */
667   for (i=0;i<deluxe_ctx->par_colors;i++) {
668     PetscSubcomm    par_subcomm;
669     const PetscInt* idxs_subproblems;
670     PetscInt        color_size;
671     PetscMPIInt     rank,active_color;
672 
673     /* get local index of i-th parallel colored problem */
674     ierr = ISGetLocalSize(is_colors[i],&color_size);CHKERRQ(ierr);
675     ierr = ISGetIndices(is_colors[i],&idxs_subproblems);CHKERRQ(ierr);
676     /* split comm for computing parallel problems for this color */
677     /* Processes not partecipating at this stage will have color = color_size */
678     /* because PetscCommDuplicate does not handle MPI_COMM_NULL */
679     active_color = color_size;
680     deluxe_ctx->par_col2sub[i] = -1;
681     for (j=0;j<n_local_parallel_problems;j++) {
682       PetscInt local_idx;
683       ierr = PetscFindInt(global_parallel[j],color_size,idxs_subproblems,&local_idx);CHKERRQ(ierr);
684       if (local_idx > -1) {
685         ierr = PetscMPIIntCast(local_idx,&active_color);CHKERRQ(ierr);
686         deluxe_ctx->par_col2sub[i] = index_parallel[j];
687         break;
688       }
689     }
690     ierr = ISRestoreIndices(is_colors[i],&idxs_subproblems);CHKERRQ(ierr);
691     ierr = PetscSubcommCreate(PetscObjectComm((PetscObject)pc),&par_subcomm);CHKERRQ(ierr);
692     ierr = PetscSubcommSetNumber(par_subcomm,color_size+1);CHKERRQ(ierr);
693     ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);CHKERRQ(ierr);
694     ierr = PetscSubcommSetTypeGeneral(par_subcomm,active_color,rank);CHKERRQ(ierr);
695     /* print debug info */
696     if (pcbddc->dbg_flag) {
697       PetscMPIInt crank,csize;
698       ierr = MPI_Comm_rank(par_subcomm->comm,&crank);CHKERRQ(ierr);
699       ierr = MPI_Comm_size(par_subcomm->comm,&csize);CHKERRQ(ierr);
700       ierr = PetscViewerASCIIPrintf(pcbddc->dbg_viewer,"Color %d: size %d, details follows.\n",i,color_size);CHKERRQ(ierr);
701       ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr);
702       ierr = PetscViewerASCIISynchronizedAllow(pcbddc->dbg_viewer,PETSC_TRUE);CHKERRQ(ierr);
703       ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"  Subdomain %d: color in subcomm %d (rank %d out of %d) (lidx %d)\n",PetscGlobalRank,par_subcomm->color,crank,csize,deluxe_ctx->par_col2sub[i]);CHKERRQ(ierr);
704       ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr);
705     }
706 
707     if (deluxe_ctx->par_col2sub[i] >= 0) {
708       PC                     pc;
709       Mat                    color_mat,color_mat_is,temp_mat;
710       ISLocalToGlobalMapping WtoNmap,l2gmap_subset;
711       IS                     is_local_numbering,isB_local,isW_local,isW;
712       PCBDDCSubSchurs        sub_schurs = deluxe_ctx->sub_schurs;
713       PetscInt               subidx,n_local_dofs,n_global_dofs;
714       PetscInt               *global_numbering,*local_numbering;
715       char                   ksp_prefix[256];
716       size_t                 len;
717 
718       /* Local index for schur complement on subset */
719       subidx = deluxe_ctx->par_col2sub[i];
720 
721       /* Parallel numbering for dofs in colored subset */
722       ierr = ISSum(sub_schurs->is_AEj_I[subidx],sub_schurs->is_AEj_B[subidx],&is_local_numbering);CHKERRQ(ierr);
723       ierr = ISGetLocalSize(is_local_numbering,&n_local_dofs);CHKERRQ(ierr);
724       ierr = ISGetIndices(is_local_numbering,(const PetscInt **)&local_numbering);CHKERRQ(ierr);
725       ierr = PCBDDCSubsetNumbering(par_subcomm->comm,pcbddc->mat_graph->l2gmap,n_local_dofs,local_numbering,PETSC_NULL,&n_global_dofs,&global_numbering);CHKERRQ(ierr);
726       ierr = ISRestoreIndices(is_local_numbering,(const PetscInt **)&local_numbering);CHKERRQ(ierr);
727 
728       /* L2Gmap from relevant dofs to local dofs */
729       ierr = ISLocalToGlobalMappingCreateIS(is_local_numbering,&WtoNmap);CHKERRQ(ierr);
730 
731       /* L2Gmap from local to global dofs */
732       ierr = ISLocalToGlobalMappingCreate(par_subcomm->comm,1,n_local_dofs,global_numbering,PETSC_COPY_VALUES,&l2gmap_subset);CHKERRQ(ierr);
733 
734       /* compute parallel matrix (extended dirichlet problem on subset) */
735       ierr = MatCreateIS(par_subcomm->comm,1,PETSC_DECIDE,PETSC_DECIDE,n_global_dofs,n_global_dofs,l2gmap_subset,&color_mat_is);CHKERRQ(ierr);
736       ierr = MatGetSubMatrix(pcbddc->local_mat,is_local_numbering,is_local_numbering,MAT_INITIAL_MATRIX,&temp_mat);CHKERRQ(ierr);
737       ierr = MatISSetLocalMat(color_mat_is,temp_mat);CHKERRQ(ierr);
738       ierr = MatDestroy(&temp_mat);CHKERRQ(ierr);
739       ierr = MatISGetMPIXAIJ(color_mat_is,MAT_INITIAL_MATRIX,&color_mat);CHKERRQ(ierr);
740       ierr = MatDestroy(&color_mat_is);CHKERRQ(ierr);
741 
742       /* work vector for (parallel) extended dirichlet problem */
743       ierr = MatCreateVecs(color_mat,&deluxe_ctx->par_vec[i],NULL);CHKERRQ(ierr);
744 
745       /* compute scatters */
746       /* deluxe_ctx->par_scctx_p[i] extension from local subset to extended dirichlet problem
747          deluxe_ctx->par_scctx_s[i] restriction from local boundary to subset -> simple copy of selected values */
748       ierr = ISGlobalToLocalMappingApplyIS(pcbddc->BtoNmap,IS_GTOLM_DROP,sub_schurs->is_AEj_B[subidx],&isB_local);CHKERRQ(ierr);
749       ierr = VecScatterCreate(pcbddc->work_scaling,isB_local,sub_schurs->work1[subidx],NULL,&deluxe_ctx->par_scctx_s[i]);CHKERRQ(ierr);
750       ierr = ISGlobalToLocalMappingApplyIS(WtoNmap,IS_GTOLM_DROP,sub_schurs->is_AEj_B[subidx],&isW_local);CHKERRQ(ierr);
751       ierr = ISLocalToGlobalMappingApplyIS(l2gmap_subset,isW_local,&isW);CHKERRQ(ierr);
752       ierr = VecScatterCreate(sub_schurs->work1[subidx],NULL,deluxe_ctx->par_vec[i],isW,&deluxe_ctx->par_scctx_p[i]);CHKERRQ(ierr);
753 
754       /* free objects no longer neeeded */
755       ierr = ISDestroy(&isW);CHKERRQ(ierr);
756       ierr = ISDestroy(&isW_local);CHKERRQ(ierr);
757       ierr = ISDestroy(&isB_local);CHKERRQ(ierr);
758       ierr = ISLocalToGlobalMappingDestroy(&WtoNmap);CHKERRQ(ierr);
759       ierr = ISLocalToGlobalMappingDestroy(&l2gmap_subset);CHKERRQ(ierr);
760       ierr = ISDestroy(&is_local_numbering);CHKERRQ(ierr);
761       ierr = PetscFree(global_numbering);CHKERRQ(ierr);
762 
763       /* KSP for extended dirichlet problem */
764       ierr = KSPCreate(par_subcomm->comm,&deluxe_ctx->par_ksp[i]);CHKERRQ(ierr);
765       ierr = KSPSetOperators(deluxe_ctx->par_ksp[i],color_mat,color_mat);CHKERRQ(ierr);
766       ierr = KSPSetTolerances(deluxe_ctx->par_ksp[i],1.e-12,1.e-12,1.e10,10000);CHKERRQ(ierr);
767       ierr = KSPSetType(deluxe_ctx->par_ksp[i],KSPPREONLY);CHKERRQ(ierr);
768       ierr = KSPGetPC(deluxe_ctx->par_ksp[i],&pc);CHKERRQ(ierr);
769       ierr = PCSetType(pc,PCREDUNDANT);CHKERRQ(ierr);
770       ierr = PetscStrlen(((PetscObject)(pcbddc->ksp_D))->prefix,&len);CHKERRQ(ierr);
771       len -= 10; /* remove "dirichlet_" */
772       ierr = PetscStrncpy(ksp_prefix,((PetscObject)(pcbddc->ksp_D))->prefix,len+1);CHKERRQ(ierr); /* PetscStrncpy puts a terminating char at the end */
773       ierr = PetscStrcat(ksp_prefix,"deluxe_par_");CHKERRQ(ierr);
774       ierr = KSPSetOptionsPrefix(deluxe_ctx->par_ksp[i],ksp_prefix);CHKERRQ(ierr);
775       ierr = KSPSetFromOptions(deluxe_ctx->par_ksp[i]);CHKERRQ(ierr);
776       ierr = KSPSetUp(deluxe_ctx->par_ksp[i]);CHKERRQ(ierr);
777       ierr = MatDestroy(&color_mat);CHKERRQ(ierr);
778     } else { /* not partecipating in color */
779       deluxe_ctx->par_ksp[i] = 0;
780       deluxe_ctx->par_vec[i] = 0;
781       deluxe_ctx->par_scctx_p[i] = 0;
782       deluxe_ctx->par_scctx_s[i] = 0;
783     }
784     ierr = PetscSubcommDestroy(&par_subcomm);CHKERRQ(ierr);
785   }
786   for (i=0;i<deluxe_ctx->par_colors;i++) {
787     ierr = ISDestroy(&is_colors[i]);CHKERRQ(ierr);
788   }
789   ierr = PetscFree(is_colors);CHKERRQ(ierr);
790 
791   if (pcbddc->dbg_flag) {
792     Vec test_vec;
793     PetscReal error;
794     PCBDDCSubSchurs sub_schurs = deluxe_ctx->sub_schurs;
795     /* test partition of unity of coloured schur complements  */
796     for (i=0;i<deluxe_ctx->par_colors;i++) {
797       PetscInt  subidx = deluxe_ctx->par_col2sub[i];
798       PetscBool error_found = PETSC_FALSE;
799       ierr = PetscViewerASCIISynchronizedAllow(pcbddc->dbg_viewer,PETSC_TRUE);CHKERRQ(ierr);
800 
801       if (deluxe_ctx->par_ksp[i]) {
802         /* create random test vec being zero on internal nodes of the extende dirichlet problem */
803         ierr = VecDuplicate(deluxe_ctx->par_vec[i],&test_vec);CHKERRQ(ierr);
804         ierr = VecSetRandom(sub_schurs->work1[subidx],PETSC_NULL);CHKERRQ(ierr);
805         ierr = VecSet(test_vec,0.0);CHKERRQ(ierr);
806         ierr = VecScatterBegin(deluxe_ctx->par_scctx_p[i],sub_schurs->work1[subidx],test_vec,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
807         ierr = VecScatterEnd(deluxe_ctx->par_scctx_p[i],sub_schurs->work1[subidx],test_vec,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
808         /* w_j */
809         ierr = VecScatterBegin(deluxe_ctx->par_scctx_p[i],test_vec,sub_schurs->work1[subidx],INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
810         ierr = VecScatterEnd(deluxe_ctx->par_scctx_p[i],test_vec,sub_schurs->work1[subidx],INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
811         /* S_j*w_j */
812         ierr = MatMult(sub_schurs->S_Ej[subidx],sub_schurs->work1[subidx],sub_schurs->work2[subidx]);CHKERRQ(ierr);
813         /* \sum_j S_j*w_j */
814         ierr = VecSet(deluxe_ctx->par_vec[i],0.0);CHKERRQ(ierr);
815         ierr = VecScatterBegin(deluxe_ctx->par_scctx_p[i],sub_schurs->work2[subidx],deluxe_ctx->par_vec[i],ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
816         ierr = VecScatterEnd(deluxe_ctx->par_scctx_p[i],sub_schurs->work2[subidx],deluxe_ctx->par_vec[i],ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
817         /* (\sum_j S_j)^(-1)(\sum_j S_j*w_j) */
818         ierr = KSPSolve(deluxe_ctx->par_ksp[i],deluxe_ctx->par_vec[i],deluxe_ctx->par_vec[i]);CHKERRQ(ierr);
819         ierr = VecScatterBegin(deluxe_ctx->par_scctx_p[i],deluxe_ctx->par_vec[i],sub_schurs->work1[subidx],INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
820         ierr = VecScatterEnd(deluxe_ctx->par_scctx_p[i],deluxe_ctx->par_vec[i],sub_schurs->work1[subidx],INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
821         ierr = VecSet(deluxe_ctx->par_vec[i],0.0);CHKERRQ(ierr);
822         ierr = VecScatterBegin(deluxe_ctx->par_scctx_p[i],sub_schurs->work1[subidx],deluxe_ctx->par_vec[i],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
823         ierr = VecScatterEnd(deluxe_ctx->par_scctx_p[i],sub_schurs->work1[subidx],deluxe_ctx->par_vec[i],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
824         /* test partition of unity */
825         ierr = VecAXPY(test_vec,-1.0,deluxe_ctx->par_vec[i]);CHKERRQ(ierr);
826         ierr = VecNorm(test_vec,NORM_INFINITY,&error);CHKERRQ(ierr);
827         if (PetscAbsReal(error) > 1.e-2) {
828           /* ierr = VecView(test_vec,0);CHKERRQ(ierr); */
829           error_found = PETSC_TRUE;
830         }
831         ierr = VecDestroy(&test_vec);CHKERRQ(ierr);
832       }
833       if (error_found) {
834         ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Error testing local schur for color %d and subdomain %d\n",i,PetscGlobalRank);CHKERRQ(ierr);
835       }
836       ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr);
837     }
838   }
839   PetscFunctionReturn(0);
840 }
841 
842 
843 #undef __FUNCT__
844 #define __FUNCT__ "PCBDDCScalingSetUp_Deluxe_Seq"
845 static PetscErrorCode PCBDDCScalingSetUp_Deluxe_Seq(PC pc,PetscInt n_local_sequential_problems,PetscInt n_sequential_problems,PetscInt global_sequential[],PetscInt local_sequential[])
846 {
847   PC_BDDC             *pcbddc=(PC_BDDC*)pc->data;
848   PCBDDCDeluxeScaling deluxe_ctx=pcbddc->deluxe_ctx;
849   PCBDDCSubSchurs     sub_schurs = deluxe_ctx->sub_schurs;
850   Mat                 global_schur_subsets,*submat_global_schur_subsets,work_mat;
851   IS                  is_to,is_from;
852   PetscScalar         *array,*fill_vals;
853   PetscInt            *all_local_idx_G,*all_local_idx_B,*all_local_idx_N,*all_permutation_G,*dummy_idx;
854   PetscInt            i,j,k,local_problem_index;
855   PetscInt            subset_size,max_subset_size,max_subset_size_red;
856   PetscInt            local_size,global_size;
857   PC                  pc_temp;
858   MatSolverPackage    solver=NULL;
859   char                ksp_prefix[256];
860   size_t              len;
861   PetscErrorCode      ierr;
862 
863   PetscFunctionBegin;
864   if (!n_sequential_problems) {
865     PetscFunctionReturn(0);
866   }
867   /* Get info on subset sizes and sum of all subsets sizes */
868   max_subset_size = 0;
869   local_size = 0;
870   for (i=0;i<n_local_sequential_problems;i++) {
871     local_problem_index = local_sequential[i];
872     ierr = ISGetLocalSize(sub_schurs->is_AEj_B[local_problem_index],&subset_size);CHKERRQ(ierr);
873     max_subset_size = PetscMax(subset_size,max_subset_size);
874     local_size += subset_size;
875   }
876 
877   /* Work arrays for local indices */
878   ierr = PetscMalloc1(local_size,&all_local_idx_B);CHKERRQ(ierr);
879   ierr = PetscMalloc1(local_size,&all_local_idx_N);CHKERRQ(ierr);
880 
881   /* Get local indices in local whole numbering and local boundary numbering */
882   local_size = 0;
883   for (i=0;i<n_local_sequential_problems;i++) {
884     PetscInt *idxs;
885     /* get info on local problem */
886     local_problem_index = local_sequential[i];
887     ierr = ISGetLocalSize(sub_schurs->is_AEj_B[local_problem_index],&subset_size);CHKERRQ(ierr);
888     ierr = ISGetIndices(sub_schurs->is_AEj_B[local_problem_index],(const PetscInt**)&idxs);CHKERRQ(ierr);
889     /* subset indices in local numbering */
890     ierr = PetscMemcpy(all_local_idx_N+local_size,idxs,subset_size*sizeof(PetscInt));CHKERRQ(ierr);
891     /* subset indices in local boundary numbering */
892     ierr = ISGlobalToLocalMappingApply(pcbddc->BtoNmap,IS_GTOLM_DROP,subset_size,idxs,&j,&all_local_idx_B[local_size]);CHKERRQ(ierr);
893     ierr = ISRestoreIndices(sub_schurs->is_AEj_B[local_problem_index],(const PetscInt**)&idxs);CHKERRQ(ierr);
894     if (j != subset_size) {
895       SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in BDDC deluxe serial %d (BtoNmap)! %d != %d\n",local_problem_index,subset_size,j);
896     }
897     local_size += subset_size;
898   }
899 
900   /* Number dofs on all subsets (parallel) and sort numbering */
901   ierr = PCBDDCSubsetNumbering(PetscObjectComm((PetscObject)pc),pcbddc->mat_graph->l2gmap,local_size,all_local_idx_N,PETSC_NULL,&global_size,&all_local_idx_G);CHKERRQ(ierr);
902   ierr = PetscMalloc1(local_size,&all_permutation_G);CHKERRQ(ierr);
903   for (i=0;i<local_size;i++) {
904     all_permutation_G[i]=i;
905   }
906   ierr = PetscSortIntWithPermutation(local_size,all_local_idx_G,all_permutation_G);CHKERRQ(ierr);
907 
908   /* Local matrix of all local Schur on subsets */
909   ierr = MatCreate(PETSC_COMM_SELF,&deluxe_ctx->seq_mat);CHKERRQ(ierr);
910   ierr = MatSetSizes(deluxe_ctx->seq_mat,PETSC_DECIDE,PETSC_DECIDE,local_size,local_size);CHKERRQ(ierr);
911   ierr = MatSetType(deluxe_ctx->seq_mat,MATAIJ);CHKERRQ(ierr);
912   ierr = MatSeqAIJSetPreallocation(deluxe_ctx->seq_mat,max_subset_size,PETSC_NULL);CHKERRQ(ierr);
913 
914   /* Global matrix of all assembled Schur on subsets */
915   ierr = MatCreate(PetscObjectComm((PetscObject)pc),&global_schur_subsets);CHKERRQ(ierr);
916   ierr = MatSetSizes(global_schur_subsets,PETSC_DECIDE,PETSC_DECIDE,global_size,global_size);CHKERRQ(ierr);
917   ierr = MatSetType(global_schur_subsets,MATAIJ);CHKERRQ(ierr);
918   ierr = MPI_Allreduce(&max_subset_size,&max_subset_size_red,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)pc));CHKERRQ(ierr);
919   ierr = MatMPIAIJSetPreallocation(global_schur_subsets,max_subset_size_red,PETSC_NULL,max_subset_size_red,PETSC_NULL);CHKERRQ(ierr);
920 
921   /* Work arrays */
922   ierr = PetscMalloc2(max_subset_size,&dummy_idx,max_subset_size*max_subset_size,&fill_vals);CHKERRQ(ierr);
923 
924   /* Loop on local problems to compute Schur complements explicitly */
925   local_size = 0;
926   for (i=0;i<n_local_sequential_problems;i++) {
927     /* get info on local problem */
928     local_problem_index = local_sequential[i];
929     ierr = ISGetLocalSize(sub_schurs->is_AEj_B[local_problem_index],&subset_size);CHKERRQ(ierr);
930     /* local Schur */
931     for (j=0;j<subset_size;j++) {
932       ierr = VecSet(sub_schurs->work1[local_problem_index],0.0);CHKERRQ(ierr);
933       ierr = VecSetValue(sub_schurs->work1[local_problem_index],j,1.0,INSERT_VALUES);CHKERRQ(ierr);
934       ierr = MatMult(sub_schurs->S_Ej[local_problem_index],sub_schurs->work1[local_problem_index],sub_schurs->work2[local_problem_index]);CHKERRQ(ierr);
935       /* store vals */
936       ierr = VecGetArray(sub_schurs->work2[local_problem_index],&array);CHKERRQ(ierr);
937       for (k=0;k<subset_size;k++) {
938         fill_vals[k*subset_size+j] = array[k];
939       }
940       ierr = VecRestoreArray(sub_schurs->work2[local_problem_index],&array);CHKERRQ(ierr);
941     }
942     for (j=0;j<subset_size;j++) {
943       dummy_idx[j]=local_size+j;
944     }
945     ierr = MatSetValues(deluxe_ctx->seq_mat,subset_size,dummy_idx,subset_size,dummy_idx,fill_vals,INSERT_VALUES);CHKERRQ(ierr);
946     ierr = MatSetValues(global_schur_subsets,subset_size,&all_local_idx_G[local_size],subset_size,&all_local_idx_G[local_size],fill_vals,ADD_VALUES);CHKERRQ(ierr);
947     local_size += subset_size;
948   }
949   ierr = MatAssemblyBegin(deluxe_ctx->seq_mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
950   ierr = MatAssemblyEnd(deluxe_ctx->seq_mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
951   ierr = MatAssemblyBegin(global_schur_subsets,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
952   ierr = MatAssemblyEnd(global_schur_subsets,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
953   ierr = PetscFree2(dummy_idx,fill_vals);CHKERRQ(ierr);
954 
955   /* Create work vectors for sequential part of deluxe */
956   ierr = MatCreateVecs(deluxe_ctx->seq_mat,&deluxe_ctx->seq_work1,&deluxe_ctx->seq_work2);CHKERRQ(ierr);
957 
958   /* Compute deluxe sequential scatter */
959   ierr = ISCreateGeneral(PETSC_COMM_SELF,local_size,all_local_idx_B,PETSC_OWN_POINTER,&is_from);CHKERRQ(ierr);
960   ierr = VecScatterCreate(pcbddc->work_scaling,is_from,deluxe_ctx->seq_work1,NULL,&deluxe_ctx->seq_scctx);CHKERRQ(ierr);
961   ierr = ISDestroy(&is_from);CHKERRQ(ierr);
962 
963   /* Get local part of (\sum_j S_Ej) */
964   for (i=0;i<local_size;i++) {
965     all_local_idx_N[i] = all_local_idx_G[all_permutation_G[i]];
966   }
967   ierr = ISCreateGeneral(PetscObjectComm((PetscObject)pc),local_size,all_local_idx_N,PETSC_OWN_POINTER,&is_to);CHKERRQ(ierr);
968   ierr = MatGetSubMatrices(global_schur_subsets,1,&is_to,&is_to,MAT_INITIAL_MATRIX,&submat_global_schur_subsets);CHKERRQ(ierr);
969   ierr = MatDestroy(&global_schur_subsets);CHKERRQ(ierr);
970   ierr = ISDestroy(&is_to);CHKERRQ(ierr);
971   for (i=0;i<local_size;i++) {
972     all_local_idx_G[all_permutation_G[i]] = i;
973   }
974   ierr = ISCreateGeneral(PETSC_COMM_SELF,local_size,all_local_idx_G,PETSC_OWN_POINTER,&is_from);CHKERRQ(ierr);
975   ierr = ISSetPermutation(is_from);CHKERRQ(ierr);
976   ierr = MatPermute(submat_global_schur_subsets[0],is_from,is_from,&work_mat);CHKERRQ(ierr);
977   ierr = MatDestroyMatrices(1,&submat_global_schur_subsets);CHKERRQ(ierr);
978   ierr = ISDestroy(&is_from);CHKERRQ(ierr);
979   ierr = PetscFree(all_permutation_G);CHKERRQ(ierr);
980 
981   /* Create KSP object for sequential part of deluxe scaling */
982   ierr = KSPCreate(PETSC_COMM_SELF,&deluxe_ctx->seq_ksp);CHKERRQ(ierr);
983   ierr = KSPSetOperators(deluxe_ctx->seq_ksp,work_mat,work_mat);CHKERRQ(ierr);
984   ierr = KSPSetType(deluxe_ctx->seq_ksp,KSPPREONLY);CHKERRQ(ierr);
985   ierr = KSPGetPC(deluxe_ctx->seq_ksp,&pc_temp);CHKERRQ(ierr);
986   ierr = PCSetType(pc_temp,PCLU);CHKERRQ(ierr);
987   ierr = KSPGetPC(pcbddc->ksp_D,&pc_temp);CHKERRQ(ierr);
988   ierr = PCFactorGetMatSolverPackage(pc_temp,(const MatSolverPackage*)&solver);CHKERRQ(ierr);
989   if (solver) {
990     PC     new_pc;
991     PCType type;
992     ierr = PCGetType(pc_temp,&type);CHKERRQ(ierr);
993     ierr = KSPGetPC(deluxe_ctx->seq_ksp,&new_pc);CHKERRQ(ierr);
994     ierr = PCSetType(new_pc,type);CHKERRQ(ierr);
995     ierr = PCFactorSetMatSolverPackage(new_pc,solver);CHKERRQ(ierr);
996   }
997   ierr = PetscStrlen(((PetscObject)(pcbddc->ksp_D))->prefix,&len);CHKERRQ(ierr);
998   len -= 10; /* remove "dirichlet_" */
999   ierr = PetscStrncpy(ksp_prefix,((PetscObject)(pcbddc->ksp_D))->prefix,len+1);CHKERRQ(ierr);
1000   ierr = PetscStrcat(ksp_prefix,"deluxe_seq_");CHKERRQ(ierr);
1001   ierr = KSPSetOptionsPrefix(deluxe_ctx->seq_ksp,ksp_prefix);CHKERRQ(ierr);
1002   ierr = KSPSetFromOptions(deluxe_ctx->seq_ksp);CHKERRQ(ierr);
1003   ierr = KSPSetUp(deluxe_ctx->seq_ksp);CHKERRQ(ierr);
1004   ierr = MatDestroy(&work_mat);CHKERRQ(ierr);
1005   PetscFunctionReturn(0);
1006 }
1007