xref: /petsc/src/ksp/pc/impls/bddc/bddcscalingbasic.c (revision 9326c5c63b46fa2f86bb5f1368bebb604ae60fb2) 
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   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingRestriction_C",NULL);CHKERRQ(ierr);
335   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingExtension_C",NULL);CHKERRQ(ierr);
336   PetscFunctionReturn(0);
337 }
338 
339 #undef __FUNCT__
340 #define __FUNCT__ "PCBDDCScalingCreate_Deluxe"
341 static PetscErrorCode PCBDDCScalingCreate_Deluxe(PC pc)
342 {
343   PC_BDDC*            pcbddc=(PC_BDDC*)pc->data;
344   PCBDDCDeluxeScaling deluxe_ctx;
345   PetscErrorCode      ierr;
346 
347   PetscFunctionBegin;
348   ierr = PetscNew(&deluxe_ctx);CHKERRQ(ierr);
349   pcbddc->deluxe_ctx = deluxe_ctx;
350   PetscFunctionReturn(0);
351 }
352 
353 #undef __FUNCT__
354 #define __FUNCT__ "PCBDDCScalingDestroy_Deluxe"
355 static PetscErrorCode PCBDDCScalingDestroy_Deluxe(PC pc)
356 {
357   PC_BDDC*            pcbddc=(PC_BDDC*)pc->data;
358   PetscErrorCode      ierr;
359 
360   PetscFunctionBegin;
361   ierr = PCBDDCScalingReset_Deluxe_Solvers(pcbddc->deluxe_ctx);CHKERRQ(ierr);
362   ierr = PetscFree(pcbddc->deluxe_ctx);CHKERRQ(ierr);
363   PetscFunctionReturn(0);
364 }
365 
366 #undef __FUNCT__
367 #define __FUNCT__ "PCBDDCScalingReset_Deluxe_Solvers"
368 static PetscErrorCode PCBDDCScalingReset_Deluxe_Solvers(PCBDDCDeluxeScaling deluxe_ctx)
369 {
370   PetscInt       i;
371   PetscErrorCode ierr;
372 
373   PetscFunctionBegin;
374   ierr = PetscFree(deluxe_ctx->idx_simple_B);CHKERRQ(ierr);
375   deluxe_ctx->n_simple = 0;
376   for (i=0;i<deluxe_ctx->seq_n;i++) {
377     ierr = VecScatterDestroy(&deluxe_ctx->seq_scctx[i]);CHKERRQ(ierr);
378     ierr = VecDestroy(&deluxe_ctx->seq_work1[i]);CHKERRQ(ierr);
379     ierr = VecDestroy(&deluxe_ctx->seq_work2[i]);CHKERRQ(ierr);
380     ierr = MatDestroy(&deluxe_ctx->seq_mat[i]);CHKERRQ(ierr);
381     ierr = MatDestroy(&deluxe_ctx->seq_mat_inv_sum[i]);CHKERRQ(ierr);
382   }
383   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);
384   ierr = PetscFree(deluxe_ctx->workspace);CHKERRQ(ierr);
385   deluxe_ctx->seq_n = 0;
386   PetscFunctionReturn(0);
387 }
388 
389 #undef __FUNCT__
390 #define __FUNCT__ "PCBDDCScalingSetUp_Deluxe"
391 static PetscErrorCode PCBDDCScalingSetUp_Deluxe(PC pc)
392 {
393   PC_IS               *pcis=(PC_IS*)pc->data;
394   PC_BDDC             *pcbddc=(PC_BDDC*)pc->data;
395   PCBDDCDeluxeScaling deluxe_ctx=pcbddc->deluxe_ctx;
396   PCBDDCSubSchurs     sub_schurs=pcbddc->sub_schurs;
397   PetscErrorCode      ierr;
398 
399   PetscFunctionBegin;
400   /* reset data structures if the topology has changed */
401   if (pcbddc->recompute_topography) {
402     ierr = PCBDDCScalingReset_Deluxe_Solvers(deluxe_ctx);CHKERRQ(ierr);
403   }
404 
405   /* Compute data structures to solve sequential problems */
406   ierr = PCBDDCScalingSetUp_Deluxe_Private(pc);CHKERRQ(ierr);
407 
408   /* diagonal scaling on interface dofs not contained in cc */
409   if (sub_schurs->is_vertices || sub_schurs->is_dir) {
410     PetscInt n_com,n_dir;
411     n_com = 0;
412     if (sub_schurs->is_vertices) {
413       ierr = ISGetLocalSize(sub_schurs->is_vertices,&n_com);CHKERRQ(ierr);
414     }
415     n_dir = 0;
416     if (sub_schurs->is_dir) {
417       ierr = ISGetLocalSize(sub_schurs->is_dir,&n_dir);CHKERRQ(ierr);
418     }
419     if (!deluxe_ctx->n_simple) {
420       deluxe_ctx->n_simple = n_dir + n_com;
421       ierr = PetscMalloc1(deluxe_ctx->n_simple,&deluxe_ctx->idx_simple_B);CHKERRQ(ierr);
422       if (sub_schurs->is_vertices) {
423         PetscInt       nmap;
424         const PetscInt *idxs;
425 
426         ierr = ISGetIndices(sub_schurs->is_vertices,&idxs);CHKERRQ(ierr);
427         ierr = ISGlobalToLocalMappingApply(pcis->BtoNmap,IS_GTOLM_DROP,n_com,idxs,&nmap,deluxe_ctx->idx_simple_B);CHKERRQ(ierr);
428         if (nmap != n_com) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error when mapping simply scaled dofs (is_vertices)! %d != %d",nmap,n_com);
429         ierr = ISRestoreIndices(sub_schurs->is_vertices,&idxs);CHKERRQ(ierr);
430       }
431       if (sub_schurs->is_dir) {
432         PetscInt       nmap;
433         const PetscInt *idxs;
434 
435         ierr = ISGetIndices(sub_schurs->is_dir,&idxs);CHKERRQ(ierr);
436         ierr = ISGlobalToLocalMappingApply(pcis->BtoNmap,IS_GTOLM_DROP,n_dir,idxs,&nmap,deluxe_ctx->idx_simple_B+n_com);CHKERRQ(ierr);
437         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);
438         ierr = ISRestoreIndices(sub_schurs->is_dir,&idxs);CHKERRQ(ierr);
439       }
440       ierr = PetscSortInt(deluxe_ctx->n_simple,deluxe_ctx->idx_simple_B);CHKERRQ(ierr);
441     } else {
442       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);
443     }
444   } else {
445     deluxe_ctx->n_simple = 0;
446     deluxe_ctx->idx_simple_B = 0;
447   }
448   PetscFunctionReturn(0);
449 }
450 
451 #undef __FUNCT__
452 #define __FUNCT__ "PCBDDCScalingSetUp_Deluxe_Private"
453 static PetscErrorCode PCBDDCScalingSetUp_Deluxe_Private(PC pc)
454 {
455   PC_BDDC                *pcbddc=(PC_BDDC*)pc->data;
456   PCBDDCDeluxeScaling    deluxe_ctx=pcbddc->deluxe_ctx;
457   PCBDDCSubSchurs        sub_schurs = pcbddc->sub_schurs;
458   PetscScalar            *matdata,*matdata2;
459   PetscInt               i,max_subset_size,cum,cum2;
460   const PetscInt         *idxs;
461   PetscBool              newsetup = PETSC_FALSE;
462   PetscErrorCode         ierr;
463 
464   PetscFunctionBegin;
465   if (!sub_schurs->n_subs) {
466     PetscFunctionReturn(0);
467   }
468 
469   /* Allocate arrays for subproblems */
470   if (!deluxe_ctx->seq_n) {
471     deluxe_ctx->seq_n = sub_schurs->n_subs;
472     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);
473     newsetup = PETSC_TRUE;
474   } else if (deluxe_ctx->seq_n != sub_schurs->n_subs) {
475     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);
476   }
477   deluxe_ctx->change_with_qr = sub_schurs->change_with_qr;
478 
479   /* Create objects for deluxe */
480   max_subset_size = 0;
481   for (i=0;i<sub_schurs->n_subs;i++) {
482     PetscInt subset_size;
483     ierr = ISGetLocalSize(sub_schurs->is_subs[i],&subset_size);CHKERRQ(ierr);
484     max_subset_size = PetscMax(subset_size,max_subset_size);
485   }
486   if (newsetup) {
487     ierr = PetscMalloc1(2*max_subset_size,&deluxe_ctx->workspace);CHKERRQ(ierr);
488   }
489   cum = cum2 = 0;
490   ierr = ISGetIndices(sub_schurs->is_Ej_all,&idxs);CHKERRQ(ierr);
491   ierr = MatSeqAIJGetArray(sub_schurs->S_Ej_all,&matdata);CHKERRQ(ierr);
492   matdata2 = NULL;
493   if (sub_schurs->sum_S_Ej_all) {
494     ierr = MatSeqAIJGetArray(sub_schurs->sum_S_Ej_all,&matdata2);CHKERRQ(ierr);
495   }
496   for (i=0;i<deluxe_ctx->seq_n;i++) {
497     PetscInt     subset_size;
498 
499     ierr = ISGetLocalSize(sub_schurs->is_subs[i],&subset_size);CHKERRQ(ierr);
500     if (newsetup) {
501       IS  sub;
502       /* work vectors */
503       ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,subset_size,deluxe_ctx->workspace,&deluxe_ctx->seq_work1[i]);CHKERRQ(ierr);
504       ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,subset_size,deluxe_ctx->workspace+subset_size,&deluxe_ctx->seq_work2[i]);CHKERRQ(ierr);
505 
506       /* scatters */
507       ierr = ISCreateGeneral(PETSC_COMM_SELF,subset_size,idxs+cum,PETSC_COPY_VALUES,&sub);CHKERRQ(ierr);
508       ierr = VecScatterCreate(pcbddc->work_scaling,sub,deluxe_ctx->seq_work1[i],NULL,&deluxe_ctx->seq_scctx[i]);CHKERRQ(ierr);
509       ierr = ISDestroy(&sub);CHKERRQ(ierr);
510 
511       /* S_E_j */
512       ierr = MatCreateSeqDense(PETSC_COMM_SELF,subset_size,subset_size,matdata+cum2,&deluxe_ctx->seq_mat[i]);CHKERRQ(ierr);
513     }
514     /* \sum_k S^k_E_j */
515     if (matdata2) {
516       ierr = MatDestroy(&deluxe_ctx->seq_mat_inv_sum[i]);CHKERRQ(ierr);
517       ierr = MatCreateSeqDense(PETSC_COMM_SELF,subset_size,subset_size,matdata2+cum2,&deluxe_ctx->seq_mat_inv_sum[i]);CHKERRQ(ierr);
518       if (sub_schurs->is_hermitian && sub_schurs->is_posdef) {
519         ierr = MatCholeskyFactor(deluxe_ctx->seq_mat_inv_sum[i],NULL,NULL);CHKERRQ(ierr);
520       } else {
521         ierr = MatLUFactor(deluxe_ctx->seq_mat_inv_sum[i],NULL,NULL,NULL);CHKERRQ(ierr);
522       }
523     }
524     cum += subset_size;
525     cum2 += subset_size*subset_size;
526   }
527   ierr = ISRestoreIndices(sub_schurs->is_Ej_all,&idxs);CHKERRQ(ierr);
528   ierr = MatSeqAIJRestoreArray(sub_schurs->S_Ej_all,&matdata);CHKERRQ(ierr);
529   if (sub_schurs->sum_S_Ej_all) {
530     ierr = MatSeqAIJRestoreArray(sub_schurs->sum_S_Ej_all,&matdata2);CHKERRQ(ierr);
531   }
532 
533   /* the change of basis is just a reference to sub_schurs->change (if any) */
534   deluxe_ctx->change = sub_schurs->change;
535   if (deluxe_ctx->change && !deluxe_ctx->change_with_qr) {
536     for (i=0;i<deluxe_ctx->seq_n;i++) {
537       if (newsetup) {
538         PC pc;
539 
540         ierr = KSPGetPC(deluxe_ctx->change[i],&pc);CHKERRQ(ierr);
541         ierr = PCSetType(pc,PCLU);CHKERRQ(ierr);
542         ierr = KSPSetFromOptions(deluxe_ctx->change[i]);CHKERRQ(ierr);
543       }
544       ierr = KSPSetUp(deluxe_ctx->change[i]);CHKERRQ(ierr);
545     }
546   }
547   PetscFunctionReturn(0);
548 }
549