xref: /petsc/src/ksp/pc/impls/bddc/bddcscalingbasic.c (revision 89d949e22d93de971d3ed7d7fdc00f23508da20b)
1 #include <bddc.h>
2 #include <bddcprivate.h>
3 
4 /* prototypes for deluxe public functions */
5 #if 0
6 extern PetscErrorCode PCBDDCScalingCreateDeluxe(PC);
7 extern PetscErrorCode PCBDDCScalingDestroyDeluxe(PC);
8 extern PetscErrorCode PCBDDCScalingSetUpDeluxe(PC);
9 #endif
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 #if 0
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   PetscScalar         *array_x,*array_D,*array;
37   PetscScalar         zero=0.0;
38   PetscInt            i;
39   PetscMPIInt         color_rank;
40   PetscErrorCode      ierr;
41 
42   /* TODO CHECK STUFF RELATED WITH FAKE WORK */
43   PetscFunctionBegin;
44   ierr = VecSet(pcbddc->work_scaling,zero);CHKERRQ(ierr); /* needed by the fake work below */
45   /* scale deluxe vertices using diagonal scaling */
46   ierr = VecGetArray(x,&array_x);CHKERRQ(ierr);
47   ierr = VecGetArray(pcis->D,&array_D);CHKERRQ(ierr);
48   ierr = VecGetArray(pcbddc->work_scaling,&array);CHKERRQ(ierr);
49   for (i=0;i<deluxe_ctx->n_simple;i++) {
50     array[deluxe_ctx->idx_simple_B[i]] = array_x[deluxe_ctx->idx_simple_B[i]]*array_D[deluxe_ctx->idx_simple_B[i]];
51   }
52   ierr = VecRestoreArray(pcbddc->work_scaling,&array);CHKERRQ(ierr);
53   ierr = VecRestoreArray(pcis->D,&array_D);CHKERRQ(ierr);
54   ierr = VecRestoreArray(x,&array_x);CHKERRQ(ierr);
55   /* sequential part : all problems and Schur applications collapsed into seq_mat at setup phase */
56   if (deluxe_ctx->seq_mat) {
57     ierr = VecScatterBegin(deluxe_ctx->seq_scctx,x,deluxe_ctx->seq_work1,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
58     ierr = VecScatterEnd(deluxe_ctx->seq_scctx,x,deluxe_ctx->seq_work1,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
59     ierr = MatMult(deluxe_ctx->seq_mat,deluxe_ctx->seq_work1,deluxe_ctx->seq_work2);CHKERRQ(ierr);
60     ierr = KSPSolve(deluxe_ctx->seq_ksp,deluxe_ctx->seq_work2,deluxe_ctx->seq_work1);CHKERRQ(ierr);
61     /* fake work due to final ADD VALUES and vertices scaling needed? TODO: check it */
62     ierr = VecScatterBegin(deluxe_ctx->seq_scctx,deluxe_ctx->seq_work1,pcbddc->work_scaling,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
63     ierr = VecScatterEnd(deluxe_ctx->seq_scctx,deluxe_ctx->seq_work1,pcbddc->work_scaling,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
64   }
65   /* parallel part */
66   for (i=0;i<deluxe_ctx->par_colors;i++) {
67     if (deluxe_ctx->par_ksp[i]) {
68       ierr = MPI_Comm_rank(deluxe_ctx->par_subcomm[i]->comm,&color_rank);CHKERRQ(ierr);
69       ierr = VecSet(deluxe_ctx->work1_B,zero);CHKERRQ(ierr);
70       ierr = VecScatterBegin(deluxe_ctx->par_scctx_s[i],x,deluxe_ctx->work1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
71       ierr = VecScatterEnd(deluxe_ctx->par_scctx_s[i],x,deluxe_ctx->work1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
72       /* apply local schur on subset S_j^-1 */
73       ierr = PCBDDCApplySchur(pc,deluxe_ctx->work1_B,deluxe_ctx->work2_B,(Vec)0,deluxe_ctx->work1_D,deluxe_ctx->work2_D);CHKERRQ(ierr);
74       /* parallel transpose solve (\sum_j S_j)^-1 */
75       ierr = VecSet(deluxe_ctx->par_vec[i],zero);CHKERRQ(ierr);
76       ierr = VecScatterBegin(deluxe_ctx->par_scctx_p[i],deluxe_ctx->work2_B,deluxe_ctx->par_vec[i],ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
77       ierr = VecScatterEnd(deluxe_ctx->par_scctx_p[i],deluxe_ctx->work2_B,deluxe_ctx->par_vec[i],ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
78       ierr = KSPSolve(deluxe_ctx->par_ksp[i],deluxe_ctx->par_vec[i],deluxe_ctx->par_vec[i]);CHKERRQ(ierr);
79       if (!color_rank) {
80         ierr = VecScatterBegin(deluxe_ctx->par_scctx_p[i],deluxe_ctx->par_vec[i],pcbddc->work_scaling,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
81         ierr = VecScatterEnd(deluxe_ctx->par_scctx_p[i],deluxe_ctx->par_vec[i],pcbddc->work_scaling,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
82       } else { /* fake work due to final ADD VALUES and vertices scaling */
83         ierr = VecScatterBegin(deluxe_ctx->par_scctx_p[i],deluxe_ctx->par_vec[i],deluxe_ctx->work1_B,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
84         ierr = VecScatterEnd(deluxe_ctx->par_scctx_p[i],deluxe_ctx->par_vec[i],deluxe_ctx->work1_B,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
85       }
86     }
87   }
88   /* put local boundary part in global vector */
89   ierr = VecSet(y,zero);CHKERRQ(ierr);
90   ierr = VecScatterBegin(pcis->global_to_B,pcbddc->work_scaling,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
91   ierr = VecScatterEnd  (pcis->global_to_B,pcbddc->work_scaling,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
92   PetscFunctionReturn(0);
93 }
94 #endif
95 
96 #undef __FUNCT__
97 #define __FUNCT__ "PCBDDCScalingExtension"
98 PetscErrorCode PCBDDCScalingExtension(PC pc, Vec local_interface_vector, Vec global_vector)
99 {
100   PC_BDDC *pcbddc=(PC_BDDC*)pc->data;
101   PetscErrorCode ierr;
102 
103   PetscFunctionBegin;
104   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
105   PetscValidHeaderSpecific(local_interface_vector,VEC_CLASSID,2);
106   PetscValidHeaderSpecific(global_vector,VEC_CLASSID,3);
107   if (local_interface_vector == pcbddc->work_scaling) {
108     SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Local vector cannot be pcbddc->work_scaling!\n");
109   }
110   ierr = PetscTryMethod(pc,"PCBDDCScalingExtension_C",(PC,Vec,Vec),(pc,local_interface_vector,global_vector));CHKERRQ(ierr);
111   PetscFunctionReturn(0);
112 }
113 
114 #undef __FUNCT__
115 #define __FUNCT__ "PCBDDCScalingRestriction_Basic"
116 static PetscErrorCode PCBDDCScalingRestriction_Basic(PC pc, Vec global_vector, Vec local_interface_vector)
117 {
118   PetscErrorCode ierr;
119   PC_IS* pcis = (PC_IS*)pc->data;
120 
121   PetscFunctionBegin;
122   ierr = VecScatterBegin(pcis->global_to_B,global_vector,local_interface_vector,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
123   ierr = VecScatterEnd(pcis->global_to_B,global_vector,local_interface_vector,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
124   /* Apply partition of unity */
125   ierr = VecPointwiseMult(local_interface_vector,pcis->D,local_interface_vector);CHKERRQ(ierr);
126   PetscFunctionReturn(0);
127 }
128 
129 #if 0
130 #undef __FUNCT__
131 #define __FUNCT__ "PCBDDCScalingRestriction_Deluxe"
132 static PetscErrorCode PCBDDCScalingRestriction_Deluxe(PC pc, Vec x, Vec y)
133 {
134   PC_IS*              pcis=(PC_IS*)pc->data;
135   PC_BDDC*            pcbddc=(PC_BDDC*)pc->data;
136   PCBDDCDeluxeScaling deluxe_ctx = pcbddc->deluxe_ctx;
137   PetscScalar         *array_y,*array_D,zero=0.0;
138   PetscInt            i;
139   PetscErrorCode      ierr;
140 
141   PetscFunctionBegin;
142   /* get local boundary part of global vector */
143   ierr = VecScatterBegin(pcis->global_to_B,x,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
144   ierr = VecScatterEnd  (pcis->global_to_B,x,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
145   /* scale vertices using diagonal scaling -> every scaling perform the same */
146   ierr = VecGetArray(y,&array_y);CHKERRQ(ierr);
147   ierr = VecGetArray(pcis->D,&array_D);CHKERRQ(ierr);
148   for (i=0;i<deluxe_ctx->n_simple;i++) {
149     array_y[deluxe_ctx->idx_simple_B[i]] *= array_D[deluxe_ctx->idx_simple_B[i]];
150   }
151   ierr = VecRestoreArray(pcis->D,&array_D);CHKERRQ(ierr);
152   ierr = VecRestoreArray(y,&array_y);CHKERRQ(ierr);
153   /* sequential part : all problems and Schur applications collapsed into seq_mat at setup phase */
154   if (deluxe_ctx->seq_mat) {
155     ierr = VecScatterBegin(deluxe_ctx->seq_scctx,y,deluxe_ctx->seq_work1,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
156     ierr = VecScatterEnd(deluxe_ctx->seq_scctx,y,deluxe_ctx->seq_work1,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
157     ierr = KSPSolveTranspose(deluxe_ctx->seq_ksp,deluxe_ctx->seq_work1,deluxe_ctx->seq_work2);CHKERRQ(ierr);
158     ierr = MatMultTranspose(deluxe_ctx->seq_mat,deluxe_ctx->seq_work2,deluxe_ctx->seq_work1);CHKERRQ(ierr);
159     ierr = VecScatterBegin(deluxe_ctx->seq_scctx,deluxe_ctx->seq_work1,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
160     ierr = VecScatterEnd(deluxe_ctx->seq_scctx,deluxe_ctx->seq_work1,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
161   }
162   /* parallel part */
163   for (i=0;i<deluxe_ctx->par_colors;i++) {
164     if (deluxe_ctx->par_ksp[i]) {
165       /* parallel solve (\sum_j S_j)^-1 */
166       ierr = VecSet(deluxe_ctx->par_vec[i],zero);CHKERRQ(ierr);
167       ierr = VecScatterBegin(deluxe_ctx->par_scctx_p[i],y,deluxe_ctx->par_vec[i],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
168       ierr = VecScatterEnd(deluxe_ctx->par_scctx_p[i],y,deluxe_ctx->par_vec[i],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
169       ierr = KSPSolveTranspose(deluxe_ctx->par_ksp[i],deluxe_ctx->par_vec[i],deluxe_ctx->par_vec[i]);CHKERRQ(ierr);
170       /* apply local schur S_j^-1 */
171       ierr = VecSet(deluxe_ctx->work1_B,zero);CHKERRQ(ierr);
172       ierr = VecScatterBegin(deluxe_ctx->par_scctx_p[i],deluxe_ctx->par_vec[i],deluxe_ctx->work1_B,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
173       ierr = VecScatterEnd(deluxe_ctx->par_scctx_p[i],deluxe_ctx->par_vec[i],deluxe_ctx->work1_B,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
174       ierr = PCBDDCApplySchurTranspose(pc,deluxe_ctx->work1_B,deluxe_ctx->work2_B,(Vec)0,deluxe_ctx->work1_D,deluxe_ctx->work2_D);CHKERRQ(ierr);
175       ierr = VecScatterBegin(deluxe_ctx->par_scctx_s[i],deluxe_ctx->work2_B,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
176       ierr = VecScatterEnd(deluxe_ctx->par_scctx_s[i],deluxe_ctx->work2_B,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
177     }
178   }
179   PetscFunctionReturn(0);
180 }
181 #endif
182 
183 #undef __FUNCT__
184 #define __FUNCT__ "PCBDDCScalingRestriction"
185 PetscErrorCode PCBDDCScalingRestriction(PC pc, Vec global_vector, Vec local_interface_vector)
186 {
187   PC_BDDC        *pcbddc=(PC_BDDC*)pc->data;
188   PetscErrorCode ierr;
189 
190   PetscFunctionBegin;
191   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
192   PetscValidHeaderSpecific(global_vector,VEC_CLASSID,2);
193   PetscValidHeaderSpecific(local_interface_vector,VEC_CLASSID,3);
194   if (local_interface_vector == pcbddc->work_scaling) {
195     SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Local vector should cannot be pcbddc->work_scaling!\n");
196   }
197   ierr = PetscTryMethod(pc,"PCBDDCScalingRestriction_C",(PC,Vec,Vec),(pc,global_vector,local_interface_vector));CHKERRQ(ierr);
198   PetscFunctionReturn(0);
199 }
200 
201 #undef __FUNCT__
202 #define __FUNCT__ "PCBDDCScalingSetUp"
203 PetscErrorCode PCBDDCScalingSetUp(PC pc)
204 {
205   PC_IS* pcis=(PC_IS*)pc->data;
206   PC_BDDC* pcbddc=(PC_BDDC*)pc->data;
207   PetscErrorCode ierr;
208 
209   PetscFunctionBegin;
210   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
211   /* create work vector for the operator */
212   ierr = VecDuplicate(pcis->vec1_B,&pcbddc->work_scaling);CHKERRQ(ierr);
213   /* rebuild pcis->D if stiffness scaling has been requested */
214   if (pcis->use_stiffness_scaling) {
215     ierr = MatGetDiagonal(pcbddc->local_mat,pcis->vec1_N);CHKERRQ(ierr);
216     ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
217     ierr = VecScatterEnd(pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
218   }
219   ierr = VecCopy(pcis->D,pcis->vec1_B);CHKERRQ(ierr);
220   ierr = VecSet(pcis->vec1_global,0.0);CHKERRQ(ierr);
221   ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
222   ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
223   ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
224   ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
225   ierr = VecPointwiseDivide(pcis->D,pcis->D,pcis->vec1_B);CHKERRQ(ierr);
226   /* now setup */
227 #if 0
228   if (pcbddc->use_deluxe_scaling) {
229     ierr = PCBDDCScalingCreateDeluxe(pc);CHKERRQ(ierr);
230     ierr = PCBDDCScalingSetUpDeluxe(pc);CHKERRQ(ierr);
231     ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingRestriction_C",PCBDDCScalingRestriction_Deluxe);CHKERRQ(ierr);
232     ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingExtension_C",PCBDDCScalingExtension_Deluxe);CHKERRQ(ierr);
233   } else {
234 #endif
235     ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingRestriction_C",PCBDDCScalingRestriction_Basic);CHKERRQ(ierr);
236     ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingExtension_C",PCBDDCScalingExtension_Basic);CHKERRQ(ierr);
237 #if 0
238   }
239 #endif
240   /* test */
241   if (pcbddc->dbg_flag) {
242     PetscViewer viewer=pcbddc->dbg_viewer;
243     PetscReal   error,gerror;
244     MPI_Comm    test_comm;
245 
246     /* extension -> from local to parallel */
247     ierr = PetscObjectGetComm((PetscObject)pc,&test_comm);CHKERRQ(ierr);
248     ierr = VecSetRandom(pcis->vec1_global,NULL);CHKERRQ(ierr);
249     ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
250     ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
251     ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,pcis->vec1_global);CHKERRQ(ierr);
252     ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_global,pcbddc->work_scaling,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
253     ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_global,pcbddc->work_scaling,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
254     ierr = VecAXPY(pcis->vec1_B,-1.0,pcbddc->work_scaling);CHKERRQ(ierr);
255     ierr = VecNorm(pcis->vec1_B,NORM_INFINITY,&error);CHKERRQ(ierr);
256     ierr = MPI_Allreduce(&error,&gerror,1,MPIU_REAL,MPI_MAX,test_comm);CHKERRQ(ierr);
257     ierr = PetscViewerASCIIPrintf(viewer,"Error scaling extension %1.14e\n",error);CHKERRQ(ierr);
258     if (PetscAbsReal(gerror)>1.e-8) {
259       ierr = VecSet(pcis->vec1_global,0.0);CHKERRQ(ierr);
260       ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
261       ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
262       ierr = VecView(pcis->vec1_global,viewer);CHKERRQ(ierr);
263     }
264     /* restriction -> from parallel to local */
265     ierr = VecSetRandom(pcis->vec1_global,NULL);CHKERRQ(ierr);
266     ierr = PCBDDCScalingRestriction(pc,pcis->vec1_global,pcis->vec1_B);CHKERRQ(ierr);
267     ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_global,pcbddc->work_scaling,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
268     ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_global,pcbddc->work_scaling,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
269     ierr = VecSet(pcis->vec1_global,0.0);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 = VecScatterBegin(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
273     ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
274     ierr = VecAXPY(pcis->vec1_B,-1.0,pcbddc->work_scaling);CHKERRQ(ierr);
275     ierr = VecNorm(pcis->vec1_B,NORM_INFINITY,&error);CHKERRQ(ierr);
276     ierr = MPI_Allreduce(&error,&gerror,1,MPIU_REAL,MPI_MAX,test_comm);CHKERRQ(ierr);
277     ierr = PetscViewerASCIIPrintf(viewer,"Error scaling restriction %1.14e\n",gerror);CHKERRQ(ierr);
278     if (PetscAbsReal(gerror)>1.e-8) {
279       ierr = VecSet(pcis->vec1_global,0.0);CHKERRQ(ierr);
280       ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
281       ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
282       ierr = VecView(pcis->vec1_global,viewer);CHKERRQ(ierr);
283     }
284   }
285   PetscFunctionReturn(0);
286 }
287 
288 #undef __FUNCT__
289 #define __FUNCT__ "PCBDDCScalingDestroy"
290 PetscErrorCode PCBDDCScalingDestroy(PC pc)
291 {
292   PC_BDDC* pcbddc=(PC_BDDC*)pc->data;
293   PetscErrorCode ierr;
294 
295   PetscFunctionBegin;
296 #if 0
297   if (pcbddc->use_deluxe_scaling) {
298     ierr = PCBDDCScalingDestroyDeluxe(pc);CHKERRQ(ierr);
299   }
300 #endif
301   ierr = VecDestroy(&pcbddc->work_scaling);CHKERRQ(ierr);
302   /* remove functions */
303   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingRestriction_C",NULL);CHKERRQ(ierr);
304   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingExtension_C",NULL);CHKERRQ(ierr);
305   PetscFunctionReturn(0);
306 }
307 
308