Actual source code: bddcscalingbasic.c
petsc-3.5.4 2015-05-23
1: #include <../src/ksp/pc/impls/bddc/bddc.h>
2: #include <../src/ksp/pc/impls/bddc/bddcprivate.h>
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
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;
20: /* Apply partition of unity */
21: VecPointwiseMult(pcbddc->work_scaling,pcis->D,local_interface_vector);
22: VecSet(global_vector,0.0);
23: VecScatterBegin(pcis->global_to_B,pcbddc->work_scaling,global_vector,ADD_VALUES,SCATTER_REVERSE);
24: VecScatterEnd(pcis->global_to_B,pcbddc->work_scaling,global_vector,ADD_VALUES,SCATTER_REVERSE);
25: return(0);
26: }
28: #if 0
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;
42: /* TODO CHECK STUFF RELATED WITH FAKE WORK */
44: VecSet(pcbddc->work_scaling,zero); /* needed by the fake work below */
45: /* scale deluxe vertices using diagonal scaling */
46: VecGetArray(x,&array_x);
47: VecGetArray(pcis->D,&array_D);
48: VecGetArray(pcbddc->work_scaling,&array);
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: VecRestoreArray(pcbddc->work_scaling,&array);
53: VecRestoreArray(pcis->D,&array_D);
54: VecRestoreArray(x,&array_x);
55: /* sequential part : all problems and Schur applications collapsed into seq_mat at setup phase */
56: if (deluxe_ctx->seq_mat) {
57: VecScatterBegin(deluxe_ctx->seq_scctx,x,deluxe_ctx->seq_work1,INSERT_VALUES,SCATTER_FORWARD);
58: VecScatterEnd(deluxe_ctx->seq_scctx,x,deluxe_ctx->seq_work1,INSERT_VALUES,SCATTER_FORWARD);
59: MatMult(deluxe_ctx->seq_mat,deluxe_ctx->seq_work1,deluxe_ctx->seq_work2);
60: KSPSolve(deluxe_ctx->seq_ksp,deluxe_ctx->seq_work2,deluxe_ctx->seq_work1);
61: /* fake work due to final ADD VALUES and vertices scaling needed? TODO: check it */
62: VecScatterBegin(deluxe_ctx->seq_scctx,deluxe_ctx->seq_work1,pcbddc->work_scaling,INSERT_VALUES,SCATTER_REVERSE);
63: VecScatterEnd(deluxe_ctx->seq_scctx,deluxe_ctx->seq_work1,pcbddc->work_scaling,INSERT_VALUES,SCATTER_REVERSE);
64: }
65: /* parallel part */
66: for (i=0;i<deluxe_ctx->par_colors;i++) {
67: if (deluxe_ctx->par_ksp[i]) {
68: MPI_Comm_rank(deluxe_ctx->par_subcomm[i]->comm,&color_rank);
69: VecSet(deluxe_ctx->work1_B,zero);
70: VecScatterBegin(deluxe_ctx->par_scctx_s[i],x,deluxe_ctx->work1_B,INSERT_VALUES,SCATTER_FORWARD);
71: VecScatterEnd(deluxe_ctx->par_scctx_s[i],x,deluxe_ctx->work1_B,INSERT_VALUES,SCATTER_FORWARD);
72: /* apply local schur on subset S_j^-1 */
73: PCBDDCApplySchur(pc,deluxe_ctx->work1_B,deluxe_ctx->work2_B,(Vec)0,deluxe_ctx->work1_D,deluxe_ctx->work2_D);
74: /* parallel transpose solve (\sum_j S_j)^-1 */
75: VecSet(deluxe_ctx->par_vec[i],zero);
76: VecScatterBegin(deluxe_ctx->par_scctx_p[i],deluxe_ctx->work2_B,deluxe_ctx->par_vec[i],ADD_VALUES,SCATTER_FORWARD);
77: VecScatterEnd(deluxe_ctx->par_scctx_p[i],deluxe_ctx->work2_B,deluxe_ctx->par_vec[i],ADD_VALUES,SCATTER_FORWARD);
78: KSPSolve(deluxe_ctx->par_ksp[i],deluxe_ctx->par_vec[i],deluxe_ctx->par_vec[i]);
79: if (!color_rank) {
80: VecScatterBegin(deluxe_ctx->par_scctx_p[i],deluxe_ctx->par_vec[i],pcbddc->work_scaling,INSERT_VALUES,SCATTER_REVERSE);
81: VecScatterEnd(deluxe_ctx->par_scctx_p[i],deluxe_ctx->par_vec[i],pcbddc->work_scaling,INSERT_VALUES,SCATTER_REVERSE);
82: } else { /* fake work due to final ADD VALUES and vertices scaling */
83: VecScatterBegin(deluxe_ctx->par_scctx_p[i],deluxe_ctx->par_vec[i],deluxe_ctx->work1_B,INSERT_VALUES,SCATTER_REVERSE);
84: VecScatterEnd(deluxe_ctx->par_scctx_p[i],deluxe_ctx->par_vec[i],deluxe_ctx->work1_B,INSERT_VALUES,SCATTER_REVERSE);
85: }
86: }
87: }
88: /* put local boundary part in global vector */
89: VecSet(y,zero);
90: VecScatterBegin(pcis->global_to_B,pcbddc->work_scaling,y,ADD_VALUES,SCATTER_REVERSE);
91: VecScatterEnd (pcis->global_to_B,pcbddc->work_scaling,y,ADD_VALUES,SCATTER_REVERSE);
92: return(0);
93: }
94: #endif
98: PetscErrorCode PCBDDCScalingExtension(PC pc, Vec local_interface_vector, Vec global_vector)
99: {
100: PC_BDDC *pcbddc=(PC_BDDC*)pc->data;
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: PetscTryMethod(pc,"PCBDDCScalingExtension_C",(PC,Vec,Vec),(pc,local_interface_vector,global_vector));
111: return(0);
112: }
116: static PetscErrorCode PCBDDCScalingRestriction_Basic(PC pc, Vec global_vector, Vec local_interface_vector)
117: {
119: PC_IS* pcis = (PC_IS*)pc->data;
122: VecScatterBegin(pcis->global_to_B,global_vector,local_interface_vector,INSERT_VALUES,SCATTER_FORWARD);
123: VecScatterEnd(pcis->global_to_B,global_vector,local_interface_vector,INSERT_VALUES,SCATTER_FORWARD);
124: /* Apply partition of unity */
125: VecPointwiseMult(local_interface_vector,pcis->D,local_interface_vector);
126: return(0);
127: }
129: #if 0
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;
142: /* get local boundary part of global vector */
143: VecScatterBegin(pcis->global_to_B,x,y,INSERT_VALUES,SCATTER_FORWARD);
144: VecScatterEnd (pcis->global_to_B,x,y,INSERT_VALUES,SCATTER_FORWARD);
145: /* scale vertices using diagonal scaling -> every scaling perform the same */
146: VecGetArray(y,&array_y);
147: VecGetArray(pcis->D,&array_D);
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: VecRestoreArray(pcis->D,&array_D);
152: VecRestoreArray(y,&array_y);
153: /* sequential part : all problems and Schur applications collapsed into seq_mat at setup phase */
154: if (deluxe_ctx->seq_mat) {
155: VecScatterBegin(deluxe_ctx->seq_scctx,y,deluxe_ctx->seq_work1,INSERT_VALUES,SCATTER_FORWARD);
156: VecScatterEnd(deluxe_ctx->seq_scctx,y,deluxe_ctx->seq_work1,INSERT_VALUES,SCATTER_FORWARD);
157: KSPSolveTranspose(deluxe_ctx->seq_ksp,deluxe_ctx->seq_work1,deluxe_ctx->seq_work2);
158: MatMultTranspose(deluxe_ctx->seq_mat,deluxe_ctx->seq_work2,deluxe_ctx->seq_work1);
159: VecScatterBegin(deluxe_ctx->seq_scctx,deluxe_ctx->seq_work1,y,INSERT_VALUES,SCATTER_REVERSE);
160: VecScatterEnd(deluxe_ctx->seq_scctx,deluxe_ctx->seq_work1,y,INSERT_VALUES,SCATTER_REVERSE);
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: VecSet(deluxe_ctx->par_vec[i],zero);
167: VecScatterBegin(deluxe_ctx->par_scctx_p[i],y,deluxe_ctx->par_vec[i],INSERT_VALUES,SCATTER_FORWARD);
168: VecScatterEnd(deluxe_ctx->par_scctx_p[i],y,deluxe_ctx->par_vec[i],INSERT_VALUES,SCATTER_FORWARD);
169: KSPSolveTranspose(deluxe_ctx->par_ksp[i],deluxe_ctx->par_vec[i],deluxe_ctx->par_vec[i]);
170: /* apply local schur S_j^-1 */
171: VecSet(deluxe_ctx->work1_B,zero);
172: VecScatterBegin(deluxe_ctx->par_scctx_p[i],deluxe_ctx->par_vec[i],deluxe_ctx->work1_B,INSERT_VALUES,SCATTER_REVERSE);
173: VecScatterEnd(deluxe_ctx->par_scctx_p[i],deluxe_ctx->par_vec[i],deluxe_ctx->work1_B,INSERT_VALUES,SCATTER_REVERSE);
174: PCBDDCApplySchurTranspose(pc,deluxe_ctx->work1_B,deluxe_ctx->work2_B,(Vec)0,deluxe_ctx->work1_D,deluxe_ctx->work2_D);
175: VecScatterBegin(deluxe_ctx->par_scctx_s[i],deluxe_ctx->work2_B,y,INSERT_VALUES,SCATTER_FORWARD);
176: VecScatterEnd(deluxe_ctx->par_scctx_s[i],deluxe_ctx->work2_B,y,INSERT_VALUES,SCATTER_FORWARD);
177: }
178: }
179: return(0);
180: }
181: #endif
185: PetscErrorCode PCBDDCScalingRestriction(PC pc, Vec global_vector, Vec local_interface_vector)
186: {
187: PC_BDDC *pcbddc=(PC_BDDC*)pc->data;
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: PetscTryMethod(pc,"PCBDDCScalingRestriction_C",(PC,Vec,Vec),(pc,global_vector,local_interface_vector));
198: return(0);
199: }
203: PetscErrorCode PCBDDCScalingSetUp(PC pc)
204: {
205: PC_IS* pcis=(PC_IS*)pc->data;
206: PC_BDDC* pcbddc=(PC_BDDC*)pc->data;
211: /* create work vector for the operator */
212: VecDuplicate(pcis->vec1_B,&pcbddc->work_scaling);
213: /* rebuild pcis->D if stiffness scaling has been requested */
214: if (pcis->use_stiffness_scaling) {
215: MatGetDiagonal(pcbddc->local_mat,pcis->vec1_N);
216: VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);
217: VecScatterEnd(pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);
218: }
219: VecCopy(pcis->D,pcis->vec1_B);
220: VecSet(pcis->vec1_global,0.0);
221: VecScatterBegin(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
222: VecScatterEnd(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
223: VecScatterBegin(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
224: VecScatterEnd(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
225: VecPointwiseDivide(pcis->D,pcis->D,pcis->vec1_B);
226: /* now setup */
227: #if 0
228: if (pcbddc->use_deluxe_scaling) {
229: PCBDDCScalingCreateDeluxe(pc);
230: PCBDDCScalingSetUpDeluxe(pc);
231: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingRestriction_C",PCBDDCScalingRestriction_Deluxe);
232: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingExtension_C",PCBDDCScalingExtension_Deluxe);
233: } else {
234: #endif
235: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingRestriction_C",PCBDDCScalingRestriction_Basic);
236: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingExtension_C",PCBDDCScalingExtension_Basic);
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;
246: /* extension -> from local to parallel */
247: PetscObjectGetComm((PetscObject)pc,&test_comm);
248: VecSetRandom(pcis->vec1_global,NULL);
249: VecScatterBegin(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
250: VecScatterEnd(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
251: PCBDDCScalingExtension(pc,pcis->vec1_B,pcis->vec1_global);
252: VecScatterBegin(pcis->global_to_B,pcis->vec1_global,pcbddc->work_scaling,INSERT_VALUES,SCATTER_FORWARD);
253: VecScatterEnd(pcis->global_to_B,pcis->vec1_global,pcbddc->work_scaling,INSERT_VALUES,SCATTER_FORWARD);
254: VecAXPY(pcis->vec1_B,-1.0,pcbddc->work_scaling);
255: VecNorm(pcis->vec1_B,NORM_INFINITY,&error);
256: MPI_Allreduce(&error,&gerror,1,MPIU_REAL,MPI_MAX,test_comm);
257: PetscViewerASCIIPrintf(viewer,"Error scaling extension %1.14e\n",error);
258: if (PetscAbsReal(gerror)>1.e-8) {
259: VecSet(pcis->vec1_global,0.0);
260: VecScatterBegin(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,INSERT_VALUES,SCATTER_REVERSE);
261: VecScatterEnd(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,INSERT_VALUES,SCATTER_REVERSE);
262: VecView(pcis->vec1_global,viewer);
263: }
264: /* restriction -> from parallel to local */
265: VecSetRandom(pcis->vec1_global,NULL);
266: PCBDDCScalingRestriction(pc,pcis->vec1_global,pcis->vec1_B);
267: VecScatterBegin(pcis->global_to_B,pcis->vec1_global,pcbddc->work_scaling,INSERT_VALUES,SCATTER_FORWARD);
268: VecScatterEnd(pcis->global_to_B,pcis->vec1_global,pcbddc->work_scaling,INSERT_VALUES,SCATTER_FORWARD);
269: VecSet(pcis->vec1_global,0.0);
270: VecScatterBegin(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
271: VecScatterEnd(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
272: VecScatterBegin(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
273: VecScatterEnd(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
274: VecAXPY(pcis->vec1_B,-1.0,pcbddc->work_scaling);
275: VecNorm(pcis->vec1_B,NORM_INFINITY,&error);
276: MPI_Allreduce(&error,&gerror,1,MPIU_REAL,MPI_MAX,test_comm);
277: PetscViewerASCIIPrintf(viewer,"Error scaling restriction %1.14e\n",gerror);
278: if (PetscAbsReal(gerror)>1.e-8) {
279: VecSet(pcis->vec1_global,0.0);
280: VecScatterBegin(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,INSERT_VALUES,SCATTER_REVERSE);
281: VecScatterEnd(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,INSERT_VALUES,SCATTER_REVERSE);
282: VecView(pcis->vec1_global,viewer);
283: }
284: }
285: return(0);
286: }
290: PetscErrorCode PCBDDCScalingDestroy(PC pc)
291: {
292: PC_BDDC* pcbddc=(PC_BDDC*)pc->data;
296: #if 0
297: if (pcbddc->use_deluxe_scaling) {
298: PCBDDCScalingDestroyDeluxe(pc);
299: }
300: #endif
301: VecDestroy(&pcbddc->work_scaling);
302: /* remove functions */
303: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingRestriction_C",NULL);
304: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingExtension_C",NULL);
305: return(0);
306: }