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