xref: /petsc/src/ksp/pc/impls/bddc/bddcscalingbasic.c (revision 35d705712e4e30cc2db45fcc137e266709f373c9) 
1 #include <../src/ksp/pc/impls/bddc/bddc.h>
2 #include <../src/ksp/pc/impls/bddc/bddcprivate.h>
3 #include <petscblaslapack.h>
4 
5 /* prototypes for deluxe functions */
6 static PetscErrorCode PCBDDCScalingCreate_Deluxe(PC);
7 static PetscErrorCode PCBDDCScalingDestroy_Deluxe(PC);
8 static PetscErrorCode PCBDDCScalingSetUp_Deluxe(PC);
9 static PetscErrorCode PCBDDCScalingSetUp_Deluxe_Private(PC);
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   PetscErrorCode      ierr;
37 
38   PetscFunctionBegin;
39   ierr = VecSet(pcbddc->work_scaling,0.0);CHKERRQ(ierr);
40   ierr = VecSet(y,0.0);CHKERRQ(ierr);
41   if (deluxe_ctx->n_simple) { /* scale deluxe vertices using diagonal scaling */
42     PetscInt          i;
43     const PetscScalar *array_x,*array_D;
44     PetscScalar       *array;
45     ierr = VecGetArrayRead(x,&array_x);CHKERRQ(ierr);
46     ierr = VecGetArrayRead(pcis->D,&array_D);CHKERRQ(ierr);
47     ierr = VecGetArray(pcbddc->work_scaling,&array);CHKERRQ(ierr);
48     for (i=0;i<deluxe_ctx->n_simple;i++) {
49       array[deluxe_ctx->idx_simple_B[i]] = array_x[deluxe_ctx->idx_simple_B[i]]*array_D[deluxe_ctx->idx_simple_B[i]];
50     }
51     ierr = VecRestoreArray(pcbddc->work_scaling,&array);CHKERRQ(ierr);
52     ierr = VecRestoreArrayRead(pcis->D,&array_D);CHKERRQ(ierr);
53     ierr = VecRestoreArrayRead(x,&array_x);CHKERRQ(ierr);
54   }
55   /* sequential part : all problems and Schur applications collapsed into a single matrix vector multiplication or a matvec and a solve */
56   if (deluxe_ctx->seq_mat) {
57     PetscInt i;
58     for (i=0;i<deluxe_ctx->seq_n;i++) {
59       if (deluxe_ctx->change) {
60         ierr = VecScatterBegin(deluxe_ctx->seq_scctx[i],x,deluxe_ctx->seq_work2[i],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
61         ierr = VecScatterEnd(deluxe_ctx->seq_scctx[i],x,deluxe_ctx->seq_work2[i],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
62         if (deluxe_ctx->change_with_qr) {
63           Mat change;
64 
65           ierr = KSPGetOperators(deluxe_ctx->change[i],&change,NULL);CHKERRQ(ierr);
66           ierr = MatMultTranspose(change,deluxe_ctx->seq_work2[i],deluxe_ctx->seq_work1[i]);CHKERRQ(ierr);
67         } else {
68           ierr = KSPSolve(deluxe_ctx->change[i],deluxe_ctx->seq_work2[i],deluxe_ctx->seq_work1[i]);CHKERRQ(ierr);
69         }
70       } else {
71         ierr = VecScatterBegin(deluxe_ctx->seq_scctx[i],x,deluxe_ctx->seq_work1[i],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
72         ierr = VecScatterEnd(deluxe_ctx->seq_scctx[i],x,deluxe_ctx->seq_work1[i],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
73       }
74       ierr = MatMultTranspose(deluxe_ctx->seq_mat[i],deluxe_ctx->seq_work1[i],deluxe_ctx->seq_work2[i]);CHKERRQ(ierr);
75       if (deluxe_ctx->seq_mat_inv_sum[i]) {
76         PetscScalar *x;
77 
78         ierr = VecGetArray(deluxe_ctx->seq_work2[i],&x);CHKERRQ(ierr);
79         ierr = VecPlaceArray(deluxe_ctx->seq_work1[i],x);CHKERRQ(ierr);
80         ierr = VecRestoreArray(deluxe_ctx->seq_work2[i],&x);CHKERRQ(ierr);
81         ierr = MatSolveTranspose(deluxe_ctx->seq_mat_inv_sum[i],deluxe_ctx->seq_work1[i],deluxe_ctx->seq_work2[i]);CHKERRQ(ierr);
82         ierr = VecResetArray(deluxe_ctx->seq_work1[i]);CHKERRQ(ierr);
83       }
84       if (deluxe_ctx->change) {
85         Mat change;
86 
87         ierr = KSPGetOperators(deluxe_ctx->change[i],&change,NULL);CHKERRQ(ierr);
88         ierr = MatMult(change,deluxe_ctx->seq_work2[i],deluxe_ctx->seq_work1[i]);CHKERRQ(ierr);
89         ierr = VecScatterBegin(deluxe_ctx->seq_scctx[i],deluxe_ctx->seq_work1[i],pcbddc->work_scaling,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
90         ierr = VecScatterEnd(deluxe_ctx->seq_scctx[i],deluxe_ctx->seq_work1[i],pcbddc->work_scaling,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
91       } else {
92         ierr = VecScatterBegin(deluxe_ctx->seq_scctx[i],deluxe_ctx->seq_work2[i],pcbddc->work_scaling,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
93         ierr = VecScatterEnd(deluxe_ctx->seq_scctx[i],deluxe_ctx->seq_work2[i],pcbddc->work_scaling,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
94       }
95     }
96   }
97   /* put local boundary part in global vector */
98   ierr = VecScatterBegin(pcis->global_to_B,pcbddc->work_scaling,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
99   ierr = VecScatterEnd(pcis->global_to_B,pcbddc->work_scaling,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
100   PetscFunctionReturn(0);
101 }
102 
103 #undef __FUNCT__
104 #define __FUNCT__ "PCBDDCScalingExtension"
105 PetscErrorCode PCBDDCScalingExtension(PC pc, Vec local_interface_vector, Vec global_vector)
106 {
107   PC_BDDC        *pcbddc=(PC_BDDC*)pc->data;
108   PetscErrorCode ierr;
109 
110   PetscFunctionBegin;
111   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
112   PetscValidHeaderSpecific(local_interface_vector,VEC_CLASSID,2);
113   PetscValidHeaderSpecific(global_vector,VEC_CLASSID,3);
114   if (local_interface_vector == pcbddc->work_scaling) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Local vector cannot be pcbddc->work_scaling!\n");
115   ierr = PetscUseMethod(pc,"PCBDDCScalingExtension_C",(PC,Vec,Vec),(pc,local_interface_vector,global_vector));CHKERRQ(ierr);
116   PetscFunctionReturn(0);
117 }
118 
119 #undef __FUNCT__
120 #define __FUNCT__ "PCBDDCScalingRestriction_Basic"
121 static PetscErrorCode PCBDDCScalingRestriction_Basic(PC pc, Vec global_vector, Vec local_interface_vector)
122 {
123   PetscErrorCode ierr;
124   PC_IS          *pcis = (PC_IS*)pc->data;
125 
126   PetscFunctionBegin;
127   ierr = VecScatterBegin(pcis->global_to_B,global_vector,local_interface_vector,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
128   ierr = VecScatterEnd(pcis->global_to_B,global_vector,local_interface_vector,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
129   /* Apply partition of unity */
130   ierr = VecPointwiseMult(local_interface_vector,pcis->D,local_interface_vector);CHKERRQ(ierr);
131   PetscFunctionReturn(0);
132 }
133 
134 #undef __FUNCT__
135 #define __FUNCT__ "PCBDDCScalingRestriction_Deluxe"
136 static PetscErrorCode PCBDDCScalingRestriction_Deluxe(PC pc, Vec x, Vec y)
137 {
138   PC_IS*              pcis=(PC_IS*)pc->data;
139   PC_BDDC*            pcbddc=(PC_BDDC*)pc->data;
140   PCBDDCDeluxeScaling deluxe_ctx = pcbddc->deluxe_ctx;
141   PetscErrorCode      ierr;
142 
143   PetscFunctionBegin;
144   /* get local boundary part of global vector */
145   ierr = VecScatterBegin(pcis->global_to_B,x,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
146   ierr = VecScatterEnd(pcis->global_to_B,x,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
147   if (deluxe_ctx->n_simple) { /* scale deluxe vertices using diagonal scaling */
148     PetscInt          i;
149     PetscScalar       *array_y;
150     const PetscScalar *array_D;
151     ierr = VecGetArray(y,&array_y);CHKERRQ(ierr);
152     ierr = VecGetArrayRead(pcis->D,&array_D);CHKERRQ(ierr);
153     for (i=0;i<deluxe_ctx->n_simple;i++) {
154       array_y[deluxe_ctx->idx_simple_B[i]] *= array_D[deluxe_ctx->idx_simple_B[i]];
155     }
156     ierr = VecRestoreArrayRead(pcis->D,&array_D);CHKERRQ(ierr);
157     ierr = VecRestoreArray(y,&array_y);CHKERRQ(ierr);
158   }
159   /* sequential part : all problems and Schur applications collapsed into a single matrix vector multiplication and ksp solution */
160   if (deluxe_ctx->seq_mat) {
161     PetscInt i;
162     for (i=0;i<deluxe_ctx->seq_n;i++) {
163       if (deluxe_ctx->change) {
164         Mat change;
165 
166         ierr = VecScatterBegin(deluxe_ctx->seq_scctx[i],y,deluxe_ctx->seq_work2[i],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
167         ierr = VecScatterEnd(deluxe_ctx->seq_scctx[i],y,deluxe_ctx->seq_work2[i],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
168         ierr = KSPGetOperators(deluxe_ctx->change[i],&change,NULL);CHKERRQ(ierr);
169         ierr = MatMultTranspose(change,deluxe_ctx->seq_work2[i],deluxe_ctx->seq_work1[i]);CHKERRQ(ierr);
170       } else {
171         ierr = VecScatterBegin(deluxe_ctx->seq_scctx[i],y,deluxe_ctx->seq_work1[i],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
172         ierr = VecScatterEnd(deluxe_ctx->seq_scctx[i],y,deluxe_ctx->seq_work1[i],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
173       }
174       if (deluxe_ctx->seq_mat_inv_sum[i]) {
175         PetscScalar *x;
176 
177         ierr = VecGetArray(deluxe_ctx->seq_work1[i],&x);CHKERRQ(ierr);
178         ierr = VecPlaceArray(deluxe_ctx->seq_work2[i],x);CHKERRQ(ierr);
179         ierr = VecRestoreArray(deluxe_ctx->seq_work1[i],&x);CHKERRQ(ierr);
180         ierr = MatSolve(deluxe_ctx->seq_mat_inv_sum[i],deluxe_ctx->seq_work1[i],deluxe_ctx->seq_work2[i]);CHKERRQ(ierr);
181         ierr = VecResetArray(deluxe_ctx->seq_work2[i]);CHKERRQ(ierr);
182       }
183       ierr = MatMult(deluxe_ctx->seq_mat[i],deluxe_ctx->seq_work1[i],deluxe_ctx->seq_work2[i]);CHKERRQ(ierr);
184       if (deluxe_ctx->change) {
185         if (deluxe_ctx->change_with_qr) {
186           Mat change;
187 
188           ierr = KSPGetOperators(deluxe_ctx->change[i],&change,NULL);CHKERRQ(ierr);
189           ierr = MatMult(change,deluxe_ctx->seq_work2[i],deluxe_ctx->seq_work1[i]);CHKERRQ(ierr);
190         } else {
191           ierr = KSPSolveTranspose(deluxe_ctx->change[i],deluxe_ctx->seq_work2[i],deluxe_ctx->seq_work1[i]);CHKERRQ(ierr);
192         }
193         ierr = VecScatterBegin(deluxe_ctx->seq_scctx[i],deluxe_ctx->seq_work1[i],y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
194         ierr = VecScatterEnd(deluxe_ctx->seq_scctx[i],deluxe_ctx->seq_work1[i],y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
195       } else {
196         ierr = VecScatterBegin(deluxe_ctx->seq_scctx[i],deluxe_ctx->seq_work2[i],y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
197         ierr = VecScatterEnd(deluxe_ctx->seq_scctx[i],deluxe_ctx->seq_work2[i],y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
198       }
199     }
200   }
201   PetscFunctionReturn(0);
202 }
203 
204 #undef __FUNCT__
205 #define __FUNCT__ "PCBDDCScalingRestriction"
206 PetscErrorCode PCBDDCScalingRestriction(PC pc, Vec global_vector, Vec local_interface_vector)
207 {
208   PC_BDDC        *pcbddc=(PC_BDDC*)pc->data;
209   PetscErrorCode ierr;
210 
211   PetscFunctionBegin;
212   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
213   PetscValidHeaderSpecific(global_vector,VEC_CLASSID,2);
214   PetscValidHeaderSpecific(local_interface_vector,VEC_CLASSID,3);
215   if (local_interface_vector == pcbddc->work_scaling) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Local vector cannot be pcbddc->work_scaling!\n");
216   ierr = PetscUseMethod(pc,"PCBDDCScalingRestriction_C",(PC,Vec,Vec),(pc,global_vector,local_interface_vector));CHKERRQ(ierr);
217   PetscFunctionReturn(0);
218 }
219 
220 #undef __FUNCT__
221 #define __FUNCT__ "PCBDDCScalingSetUp"
222 PetscErrorCode PCBDDCScalingSetUp(PC pc)
223 {
224   PC_IS*         pcis=(PC_IS*)pc->data;
225   PC_BDDC*       pcbddc=(PC_BDDC*)pc->data;
226   PetscErrorCode ierr;
227 
228   PetscFunctionBegin;
229   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
230   /* create work vector for the operator */
231   ierr = VecDestroy(&pcbddc->work_scaling);CHKERRQ(ierr);
232   ierr = VecDuplicate(pcis->vec1_B,&pcbddc->work_scaling);CHKERRQ(ierr);
233   /* always rebuild pcis->D */
234   if (pcis->use_stiffness_scaling) {
235     ierr = MatGetDiagonal(pcbddc->local_mat,pcis->vec1_N);CHKERRQ(ierr);
236     ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
237     ierr = VecScatterEnd(pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
238   }
239   ierr = VecCopy(pcis->D,pcis->vec1_B);CHKERRQ(ierr);
240   ierr = VecSet(pcis->vec1_global,0.0);CHKERRQ(ierr);
241   ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
242   ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
243   ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
244   ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
245   ierr = VecPointwiseDivide(pcis->D,pcis->D,pcis->vec1_B);CHKERRQ(ierr);
246   /* now setup */
247   if (pcbddc->use_deluxe_scaling) {
248     if (!pcbddc->deluxe_ctx) {
249       ierr = PCBDDCScalingCreate_Deluxe(pc);CHKERRQ(ierr);
250     }
251     ierr = PCBDDCScalingSetUp_Deluxe(pc);CHKERRQ(ierr);
252     ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingRestriction_C",PCBDDCScalingRestriction_Deluxe);CHKERRQ(ierr);
253     ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingExtension_C",PCBDDCScalingExtension_Deluxe);CHKERRQ(ierr);
254   } else {
255     ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingRestriction_C",PCBDDCScalingRestriction_Basic);CHKERRQ(ierr);
256     ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingExtension_C",PCBDDCScalingExtension_Basic);CHKERRQ(ierr);
257   }
258 
259   /* test */
260   if (pcbddc->dbg_flag) {
261     Mat         B0_B = NULL;
262     Vec         B0_Bv = NULL, B0_Bv2 = NULL;
263     Vec         vec2_global;
264     PetscViewer viewer = pcbddc->dbg_viewer;
265     PetscReal   error;
266 
267     /* extension -> from local to parallel */
268     ierr = VecSet(pcis->vec1_global,0.0);CHKERRQ(ierr);
269     ierr = VecSetRandom(pcis->vec1_B,NULL);CHKERRQ(ierr);
270     ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
271     ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
272     ierr = VecDuplicate(pcis->vec1_global,&vec2_global);CHKERRQ(ierr);
273     ierr = VecCopy(pcis->vec1_global,vec2_global);CHKERRQ(ierr);
274     ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
275     ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
276     if (pcbddc->benign_n) {
277       IS is_dummy;
278 
279       ierr = ISCreateStride(PETSC_COMM_SELF,pcbddc->benign_n,0,1,&is_dummy);CHKERRQ(ierr);
280       ierr = MatGetSubMatrix(pcbddc->benign_B0,is_dummy,pcis->is_B_local,MAT_INITIAL_MATRIX,&B0_B);CHKERRQ(ierr);
281       ierr = ISDestroy(&is_dummy);CHKERRQ(ierr);
282       ierr = MatCreateVecs(B0_B,NULL,&B0_Bv);CHKERRQ(ierr);
283       ierr = VecDuplicate(B0_Bv,&B0_Bv2);CHKERRQ(ierr);
284       ierr = MatMult(B0_B,pcis->vec1_B,B0_Bv);CHKERRQ(ierr);
285     }
286     ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,pcis->vec1_global);CHKERRQ(ierr);
287     if (pcbddc->benign_saddle_point) {
288       PetscReal errorl = 0.;
289       ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
290       ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
291       if (pcbddc->benign_n) {
292         ierr = MatMult(B0_B,pcis->vec1_B,B0_Bv2);CHKERRQ(ierr);
293         ierr = VecAXPY(B0_Bv,-1.0,B0_Bv2);CHKERRQ(ierr);
294         ierr = VecNorm(B0_Bv,NORM_INFINITY,&errorl);CHKERRQ(ierr);
295       }
296       ierr = MPI_Allreduce(&errorl,&error,1,MPIU_REAL,MPI_SUM,PetscObjectComm((PetscObject)pc));CHKERRQ(ierr);
297       ierr = PetscViewerASCIIPrintf(viewer,"Error benign extension %1.14e\n",error);CHKERRQ(ierr);
298     }
299     ierr = VecAXPY(pcis->vec1_global,-1.0,vec2_global);CHKERRQ(ierr);
300     ierr = VecNorm(pcis->vec1_global,NORM_INFINITY,&error);CHKERRQ(ierr);
301     ierr = PetscViewerASCIIPrintf(viewer,"Error scaling extension %1.14e\n",error);CHKERRQ(ierr);
302     ierr = VecDestroy(&vec2_global);CHKERRQ(ierr);
303 
304     /* restriction -> from parallel to local */
305     ierr = VecSet(pcis->vec1_global,0.0);CHKERRQ(ierr);
306     ierr = VecSetRandom(pcis->vec1_B,NULL);CHKERRQ(ierr);
307     ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
308     ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
309     ierr = PCBDDCScalingRestriction(pc,pcis->vec1_global,pcis->vec1_B);CHKERRQ(ierr);
310     ierr = VecScale(pcis->vec1_B,-1.0);CHKERRQ(ierr);
311     ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
312     ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
313     ierr = VecNorm(pcis->vec1_global,NORM_INFINITY,&error);CHKERRQ(ierr);
314     ierr = PetscViewerASCIIPrintf(viewer,"Error scaling restriction %1.14e\n",error);CHKERRQ(ierr);
315     ierr = MatDestroy(&B0_B);CHKERRQ(ierr);
316     ierr = VecDestroy(&B0_Bv);CHKERRQ(ierr);
317     ierr = VecDestroy(&B0_Bv2);CHKERRQ(ierr);
318   }
319   PetscFunctionReturn(0);
320 }
321 
322 #undef __FUNCT__
323 #define __FUNCT__ "PCBDDCScalingDestroy"
324 PetscErrorCode PCBDDCScalingDestroy(PC pc)
325 {
326   PC_BDDC* pcbddc=(PC_BDDC*)pc->data;
327   PetscErrorCode ierr;
328 
329   PetscFunctionBegin;
330   if (pcbddc->deluxe_ctx) {
331     ierr = PCBDDCScalingDestroy_Deluxe(pc);CHKERRQ(ierr);
332   }
333   ierr = VecDestroy(&pcbddc->work_scaling);CHKERRQ(ierr);
334   /* remove functions */
335   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingRestriction_C",NULL);CHKERRQ(ierr);
336   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingExtension_C",NULL);CHKERRQ(ierr);
337   PetscFunctionReturn(0);
338 }
339 
340 #undef __FUNCT__
341 #define __FUNCT__ "PCBDDCScalingCreate_Deluxe"
342 static PetscErrorCode PCBDDCScalingCreate_Deluxe(PC pc)
343 {
344   PC_BDDC*            pcbddc=(PC_BDDC*)pc->data;
345   PCBDDCDeluxeScaling deluxe_ctx;
346   PetscErrorCode      ierr;
347 
348   PetscFunctionBegin;
349   ierr = PetscNew(&deluxe_ctx);CHKERRQ(ierr);
350   pcbddc->deluxe_ctx = deluxe_ctx;
351   PetscFunctionReturn(0);
352 }
353 
354 #undef __FUNCT__
355 #define __FUNCT__ "PCBDDCScalingDestroy_Deluxe"
356 static PetscErrorCode PCBDDCScalingDestroy_Deluxe(PC pc)
357 {
358   PC_BDDC*            pcbddc=(PC_BDDC*)pc->data;
359   PetscErrorCode      ierr;
360 
361   PetscFunctionBegin;
362   ierr = PCBDDCScalingReset_Deluxe_Solvers(pcbddc->deluxe_ctx);CHKERRQ(ierr);
363   ierr = PetscFree(pcbddc->deluxe_ctx);CHKERRQ(ierr);
364   PetscFunctionReturn(0);
365 }
366 
367 #undef __FUNCT__
368 #define __FUNCT__ "PCBDDCScalingReset_Deluxe_Solvers"
369 static PetscErrorCode PCBDDCScalingReset_Deluxe_Solvers(PCBDDCDeluxeScaling deluxe_ctx)
370 {
371   PetscInt       i;
372   PetscErrorCode ierr;
373 
374   PetscFunctionBegin;
375   ierr = PetscFree(deluxe_ctx->idx_simple_B);CHKERRQ(ierr);
376   deluxe_ctx->n_simple = 0;
377   for (i=0;i<deluxe_ctx->seq_n;i++) {
378     ierr = VecScatterDestroy(&deluxe_ctx->seq_scctx[i]);CHKERRQ(ierr);
379     ierr = VecDestroy(&deluxe_ctx->seq_work1[i]);CHKERRQ(ierr);
380     ierr = VecDestroy(&deluxe_ctx->seq_work2[i]);CHKERRQ(ierr);
381     ierr = MatDestroy(&deluxe_ctx->seq_mat[i]);CHKERRQ(ierr);
382     ierr = MatDestroy(&deluxe_ctx->seq_mat_inv_sum[i]);CHKERRQ(ierr);
383   }
384   ierr = PetscFree5(deluxe_ctx->seq_scctx,deluxe_ctx->seq_work1,deluxe_ctx->seq_work2,deluxe_ctx->seq_mat,deluxe_ctx->seq_mat_inv_sum);CHKERRQ(ierr);
385   ierr = PetscFree(deluxe_ctx->workspace);CHKERRQ(ierr);
386   deluxe_ctx->seq_n = 0;
387   PetscFunctionReturn(0);
388 }
389 
390 #undef __FUNCT__
391 #define __FUNCT__ "PCBDDCScalingSetUp_Deluxe"
392 static PetscErrorCode PCBDDCScalingSetUp_Deluxe(PC pc)
393 {
394   PC_IS               *pcis=(PC_IS*)pc->data;
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   /* reset data structures if the topology has changed */
402   if (pcbddc->recompute_topography) {
403     ierr = PCBDDCScalingReset_Deluxe_Solvers(deluxe_ctx);CHKERRQ(ierr);
404   }
405 
406   /* Compute data structures to solve sequential problems */
407   ierr = PCBDDCScalingSetUp_Deluxe_Private(pc);CHKERRQ(ierr);
408 
409   /* diagonal scaling on interface dofs not contained in cc */
410   if (sub_schurs->is_vertices || sub_schurs->is_dir) {
411     PetscInt n_com,n_dir;
412     n_com = 0;
413     if (sub_schurs->is_vertices) {
414       ierr = ISGetLocalSize(sub_schurs->is_vertices,&n_com);CHKERRQ(ierr);
415     }
416     n_dir = 0;
417     if (sub_schurs->is_dir) {
418       ierr = ISGetLocalSize(sub_schurs->is_dir,&n_dir);CHKERRQ(ierr);
419     }
420     if (!deluxe_ctx->n_simple) {
421       deluxe_ctx->n_simple = n_dir + n_com;
422       ierr = PetscMalloc1(deluxe_ctx->n_simple,&deluxe_ctx->idx_simple_B);CHKERRQ(ierr);
423       if (sub_schurs->is_vertices) {
424         PetscInt       nmap;
425         const PetscInt *idxs;
426 
427         ierr = ISGetIndices(sub_schurs->is_vertices,&idxs);CHKERRQ(ierr);
428         ierr = ISGlobalToLocalMappingApply(pcis->BtoNmap,IS_GTOLM_DROP,n_com,idxs,&nmap,deluxe_ctx->idx_simple_B);CHKERRQ(ierr);
429         if (nmap != n_com) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error when mapping simply scaled dofs (is_vertices)! %d != %d",nmap,n_com);
430         ierr = ISRestoreIndices(sub_schurs->is_vertices,&idxs);CHKERRQ(ierr);
431       }
432       if (sub_schurs->is_dir) {
433         PetscInt       nmap;
434         const PetscInt *idxs;
435 
436         ierr = ISGetIndices(sub_schurs->is_dir,&idxs);CHKERRQ(ierr);
437         ierr = ISGlobalToLocalMappingApply(pcis->BtoNmap,IS_GTOLM_DROP,n_dir,idxs,&nmap,deluxe_ctx->idx_simple_B+n_com);CHKERRQ(ierr);
438         if (nmap != n_dir) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error when mapping simply scaled dofs (sub_schurs->is_dir)! %d != %d",nmap,n_dir);
439         ierr = ISRestoreIndices(sub_schurs->is_dir,&idxs);CHKERRQ(ierr);
440       }
441       ierr = PetscSortInt(deluxe_ctx->n_simple,deluxe_ctx->idx_simple_B);CHKERRQ(ierr);
442     } else {
443       if (deluxe_ctx->n_simple != n_dir + n_com) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Number of simply scaled dofs %d is different from the previous one computed %d",n_dir + n_com,deluxe_ctx->n_simple);
444     }
445   } else {
446     deluxe_ctx->n_simple = 0;
447     deluxe_ctx->idx_simple_B = 0;
448   }
449   PetscFunctionReturn(0);
450 }
451 
452 #undef __FUNCT__
453 #define __FUNCT__ "PCBDDCScalingSetUp_Deluxe_Private"
454 static PetscErrorCode PCBDDCScalingSetUp_Deluxe_Private(PC pc)
455 {
456   PC_BDDC                *pcbddc=(PC_BDDC*)pc->data;
457   PCBDDCDeluxeScaling    deluxe_ctx=pcbddc->deluxe_ctx;
458   PCBDDCSubSchurs        sub_schurs = pcbddc->sub_schurs;
459   PetscScalar            *matdata,*matdata2;
460   PetscInt               i,max_subset_size,cum,cum2;
461   const PetscInt         *idxs;
462   PetscBool              newsetup = PETSC_FALSE;
463   PetscErrorCode         ierr;
464 
465   PetscFunctionBegin;
466   if (!sub_schurs->n_subs) {
467     PetscFunctionReturn(0);
468   }
469 
470   /* Allocate arrays for subproblems */
471   if (!deluxe_ctx->seq_n) {
472     deluxe_ctx->seq_n = sub_schurs->n_subs;
473     ierr = PetscCalloc5(deluxe_ctx->seq_n,&deluxe_ctx->seq_scctx,deluxe_ctx->seq_n,&deluxe_ctx->seq_work1,deluxe_ctx->seq_n,&deluxe_ctx->seq_work2,deluxe_ctx->seq_n,&deluxe_ctx->seq_mat,deluxe_ctx->seq_n,&deluxe_ctx->seq_mat_inv_sum);CHKERRQ(ierr);
474     newsetup = PETSC_TRUE;
475   } else if (deluxe_ctx->seq_n != sub_schurs->n_subs) {
476     SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Number of deluxe subproblems %d is different from the sub_schurs %d",deluxe_ctx->seq_n,sub_schurs->n_subs);
477   }
478   deluxe_ctx->change_with_qr = sub_schurs->change_with_qr;
479 
480   /* Create objects for deluxe */
481   max_subset_size = 0;
482   for (i=0;i<sub_schurs->n_subs;i++) {
483     PetscInt subset_size;
484     ierr = ISGetLocalSize(sub_schurs->is_subs[i],&subset_size);CHKERRQ(ierr);
485     max_subset_size = PetscMax(subset_size,max_subset_size);
486   }
487   if (newsetup) {
488     ierr = PetscMalloc1(2*max_subset_size,&deluxe_ctx->workspace);CHKERRQ(ierr);
489   }
490   cum = cum2 = 0;
491   ierr = ISGetIndices(sub_schurs->is_Ej_all,&idxs);CHKERRQ(ierr);
492   ierr = MatSeqAIJGetArray(sub_schurs->S_Ej_all,&matdata);CHKERRQ(ierr);
493   matdata2 = NULL;
494   if (sub_schurs->sum_S_Ej_all) {
495     ierr = MatSeqAIJGetArray(sub_schurs->sum_S_Ej_all,&matdata2);CHKERRQ(ierr);
496   }
497   for (i=0;i<deluxe_ctx->seq_n;i++) {
498     PetscInt     subset_size;
499 
500     ierr = ISGetLocalSize(sub_schurs->is_subs[i],&subset_size);CHKERRQ(ierr);
501     if (newsetup) {
502       IS  sub;
503       /* work vectors */
504       ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,subset_size,deluxe_ctx->workspace,&deluxe_ctx->seq_work1[i]);CHKERRQ(ierr);
505       ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,subset_size,deluxe_ctx->workspace+subset_size,&deluxe_ctx->seq_work2[i]);CHKERRQ(ierr);
506 
507       /* scatters */
508       ierr = ISCreateGeneral(PETSC_COMM_SELF,subset_size,idxs+cum,PETSC_COPY_VALUES,&sub);CHKERRQ(ierr);
509       ierr = VecScatterCreate(pcbddc->work_scaling,sub,deluxe_ctx->seq_work1[i],NULL,&deluxe_ctx->seq_scctx[i]);CHKERRQ(ierr);
510       ierr = ISDestroy(&sub);CHKERRQ(ierr);
511 
512       /* S_E_j */
513       ierr = MatCreateSeqDense(PETSC_COMM_SELF,subset_size,subset_size,matdata+cum2,&deluxe_ctx->seq_mat[i]);CHKERRQ(ierr);
514     }
515     /* \sum_k S^k_E_j */
516     if (matdata2) {
517       ierr = MatDestroy(&deluxe_ctx->seq_mat_inv_sum[i]);CHKERRQ(ierr);
518       ierr = MatCreateSeqDense(PETSC_COMM_SELF,subset_size,subset_size,matdata2+cum2,&deluxe_ctx->seq_mat_inv_sum[i]);CHKERRQ(ierr);
519       if (sub_schurs->is_hermitian && sub_schurs->is_posdef) {
520         ierr = MatCholeskyFactor(deluxe_ctx->seq_mat_inv_sum[i],NULL,NULL);CHKERRQ(ierr);
521       } else {
522         ierr = MatLUFactor(deluxe_ctx->seq_mat_inv_sum[i],NULL,NULL,NULL);CHKERRQ(ierr);
523       }
524     }
525     cum += subset_size;
526     cum2 += subset_size*subset_size;
527   }
528   ierr = ISRestoreIndices(sub_schurs->is_Ej_all,&idxs);CHKERRQ(ierr);
529   ierr = MatSeqAIJRestoreArray(sub_schurs->S_Ej_all,&matdata);CHKERRQ(ierr);
530   if (sub_schurs->sum_S_Ej_all) {
531     ierr = MatSeqAIJRestoreArray(sub_schurs->sum_S_Ej_all,&matdata2);CHKERRQ(ierr);
532   }
533 
534   /* the change of basis is just a reference to sub_schurs->change (if any) */
535   deluxe_ctx->change = sub_schurs->change;
536   if (deluxe_ctx->change && !deluxe_ctx->change_with_qr) {
537     for (i=0;i<deluxe_ctx->seq_n;i++) {
538       if (newsetup) {
539         PC pc;
540 
541         ierr = KSPGetPC(deluxe_ctx->change[i],&pc);CHKERRQ(ierr);
542         ierr = PCSetType(pc,PCLU);CHKERRQ(ierr);
543         ierr = KSPSetFromOptions(deluxe_ctx->change[i]);CHKERRQ(ierr);
544       }
545       ierr = KSPSetUp(deluxe_ctx->change[i]);CHKERRQ(ierr);
546     }
547   }
548   PetscFunctionReturn(0);
549 }
550