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