Actual source code: bddcscalingbasic.c
petsc-3.6.4 2016-04-12
1: #include <../src/ksp/pc/impls/bddc/bddc.h>
2: #include <../src/ksp/pc/impls/bddc/bddcprivate.h>
4: /* prototypes for deluxe functions */
5: static PetscErrorCode PCBDDCScalingCreate_Deluxe(PC);
6: static PetscErrorCode PCBDDCScalingDestroy_Deluxe(PC);
7: static PetscErrorCode PCBDDCScalingSetUp_Deluxe(PC);
8: static PetscErrorCode PCBDDCScalingSetUp_Deluxe_Private(PC);
9: static PetscErrorCode PCBDDCScalingReset_Deluxe_Solvers(PCBDDCDeluxeScaling);
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: }
30: static PetscErrorCode PCBDDCScalingExtension_Deluxe(PC pc, Vec x, Vec y)
31: {
32: PC_IS* pcis=(PC_IS*)pc->data;
33: PC_BDDC* pcbddc=(PC_BDDC*)pc->data;
34: PCBDDCDeluxeScaling deluxe_ctx = pcbddc->deluxe_ctx;
35: PetscErrorCode ierr;
38: VecSet(pcbddc->work_scaling,0.0);
39: VecSet(y,0.0);
40: if (deluxe_ctx->n_simple) { /* scale deluxe vertices using diagonal scaling */
41: PetscInt i;
42: const PetscScalar *array_x,*array_D;
43: PetscScalar *array;
44: VecGetArrayRead(x,&array_x);
45: VecGetArrayRead(pcis->D,&array_D);
46: VecGetArray(pcbddc->work_scaling,&array);
47: for (i=0;i<deluxe_ctx->n_simple;i++) {
48: array[deluxe_ctx->idx_simple_B[i]] = array_x[deluxe_ctx->idx_simple_B[i]]*array_D[deluxe_ctx->idx_simple_B[i]];
49: }
50: VecRestoreArray(pcbddc->work_scaling,&array);
51: VecRestoreArrayRead(pcis->D,&array_D);
52: VecRestoreArrayRead(x,&array_x);
53: }
54: /* sequential part : all problems and Schur applications collapsed into a single matrix vector multiplication or a matvec and a solve */
55: if (deluxe_ctx->seq_mat) {
56: VecScatterBegin(deluxe_ctx->seq_scctx,x,deluxe_ctx->seq_work1,INSERT_VALUES,SCATTER_FORWARD);
57: VecScatterEnd(deluxe_ctx->seq_scctx,x,deluxe_ctx->seq_work1,INSERT_VALUES,SCATTER_FORWARD);
58: MatMultTranspose(deluxe_ctx->seq_mat,deluxe_ctx->seq_work1,deluxe_ctx->seq_work2);
59: if (deluxe_ctx->seq_ksp) {
60: KSPSolveTranspose(deluxe_ctx->seq_ksp,deluxe_ctx->seq_work2,deluxe_ctx->seq_work2);
61: }
62: VecScatterBegin(deluxe_ctx->seq_scctx,deluxe_ctx->seq_work2,pcbddc->work_scaling,INSERT_VALUES,SCATTER_REVERSE);
63: VecScatterEnd(deluxe_ctx->seq_scctx,deluxe_ctx->seq_work2,pcbddc->work_scaling,INSERT_VALUES,SCATTER_REVERSE);
64: }
65: /* put local boundary part in global vector */
66: VecScatterBegin(pcis->global_to_B,pcbddc->work_scaling,y,ADD_VALUES,SCATTER_REVERSE);
67: VecScatterEnd(pcis->global_to_B,pcbddc->work_scaling,y,ADD_VALUES,SCATTER_REVERSE);
68: return(0);
69: }
73: PetscErrorCode PCBDDCScalingExtension(PC pc, Vec local_interface_vector, Vec global_vector)
74: {
75: PC_BDDC *pcbddc=(PC_BDDC*)pc->data;
82: if (local_interface_vector == pcbddc->work_scaling) {
83: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Local vector cannot be pcbddc->work_scaling!\n");
84: }
85: PetscTryMethod(pc,"PCBDDCScalingExtension_C",(PC,Vec,Vec),(pc,local_interface_vector,global_vector));
86: return(0);
87: }
91: static PetscErrorCode PCBDDCScalingRestriction_Basic(PC pc, Vec global_vector, Vec local_interface_vector)
92: {
94: PC_IS* pcis = (PC_IS*)pc->data;
97: VecScatterBegin(pcis->global_to_B,global_vector,local_interface_vector,INSERT_VALUES,SCATTER_FORWARD);
98: VecScatterEnd(pcis->global_to_B,global_vector,local_interface_vector,INSERT_VALUES,SCATTER_FORWARD);
99: /* Apply partition of unity */
100: VecPointwiseMult(local_interface_vector,pcis->D,local_interface_vector);
101: return(0);
102: }
106: static PetscErrorCode PCBDDCScalingRestriction_Deluxe(PC pc, Vec x, Vec y)
107: {
108: PC_IS* pcis=(PC_IS*)pc->data;
109: PC_BDDC* pcbddc=(PC_BDDC*)pc->data;
110: PCBDDCDeluxeScaling deluxe_ctx = pcbddc->deluxe_ctx;
111: PetscErrorCode ierr;
114: /* get local boundary part of global vector */
115: VecScatterBegin(pcis->global_to_B,x,y,INSERT_VALUES,SCATTER_FORWARD);
116: VecScatterEnd(pcis->global_to_B,x,y,INSERT_VALUES,SCATTER_FORWARD);
117: if (deluxe_ctx->n_simple) { /* scale deluxe vertices using diagonal scaling */
118: PetscInt i;
119: PetscScalar *array_y;
120: const PetscScalar *array_D;
121: VecGetArray(y,&array_y);
122: VecGetArrayRead(pcis->D,&array_D);
123: for (i=0;i<deluxe_ctx->n_simple;i++) {
124: array_y[deluxe_ctx->idx_simple_B[i]] *= array_D[deluxe_ctx->idx_simple_B[i]];
125: }
126: VecRestoreArrayRead(pcis->D,&array_D);
127: VecRestoreArray(y,&array_y);
128: }
129: /* sequential part : all problems and Schur applications collapsed into a single matrix vector multiplication and ksp solution */
130: if (deluxe_ctx->seq_mat) {
131: VecScatterBegin(deluxe_ctx->seq_scctx,y,deluxe_ctx->seq_work1,INSERT_VALUES,SCATTER_FORWARD);
132: VecScatterEnd(deluxe_ctx->seq_scctx,y,deluxe_ctx->seq_work1,INSERT_VALUES,SCATTER_FORWARD);
133: if (deluxe_ctx->seq_ksp) {
134: KSPSolve(deluxe_ctx->seq_ksp,deluxe_ctx->seq_work1,deluxe_ctx->seq_work1);
135: }
136: MatMult(deluxe_ctx->seq_mat,deluxe_ctx->seq_work1,deluxe_ctx->seq_work2);
137: VecScatterBegin(deluxe_ctx->seq_scctx,deluxe_ctx->seq_work2,y,INSERT_VALUES,SCATTER_REVERSE);
138: VecScatterEnd(deluxe_ctx->seq_scctx,deluxe_ctx->seq_work2,y,INSERT_VALUES,SCATTER_REVERSE);
139: }
140: return(0);
141: }
145: PetscErrorCode PCBDDCScalingRestriction(PC pc, Vec global_vector, Vec local_interface_vector)
146: {
147: PC_BDDC *pcbddc=(PC_BDDC*)pc->data;
154: if (local_interface_vector == pcbddc->work_scaling) {
155: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Local vector cannot be pcbddc->work_scaling!\n");
156: }
157: PetscTryMethod(pc,"PCBDDCScalingRestriction_C",(PC,Vec,Vec),(pc,global_vector,local_interface_vector));
158: return(0);
159: }
163: PetscErrorCode PCBDDCScalingSetUp(PC pc)
164: {
165: PC_IS* pcis=(PC_IS*)pc->data;
166: PC_BDDC* pcbddc=(PC_BDDC*)pc->data;
171: /* create work vector for the operator */
172: VecDestroy(&pcbddc->work_scaling);
173: VecDuplicate(pcis->vec1_B,&pcbddc->work_scaling);
174: /* always rebuild pcis->D */
175: if (pcis->use_stiffness_scaling) {
176: MatGetDiagonal(pcbddc->local_mat,pcis->vec1_N);
177: VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);
178: VecScatterEnd(pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);
179: }
180: VecCopy(pcis->D,pcis->vec1_B);
181: VecSet(pcis->vec1_global,0.0);
182: VecScatterBegin(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
183: VecScatterEnd(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
184: VecScatterBegin(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
185: VecScatterEnd(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
186: VecPointwiseDivide(pcis->D,pcis->D,pcis->vec1_B);
187: /* now setup */
188: if (pcbddc->use_deluxe_scaling) {
189: if (!pcbddc->deluxe_ctx) {
190: PCBDDCScalingCreate_Deluxe(pc);
191: }
192: PCBDDCScalingSetUp_Deluxe(pc);
193: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingRestriction_C",PCBDDCScalingRestriction_Deluxe);
194: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingExtension_C",PCBDDCScalingExtension_Deluxe);
195: } else {
196: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingRestriction_C",PCBDDCScalingRestriction_Basic);
197: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingExtension_C",PCBDDCScalingExtension_Basic);
198: }
200: /* test */
201: if (pcbddc->dbg_flag) {
202: Vec vec2_global;
203: PetscViewer viewer=pcbddc->dbg_viewer;
204: PetscReal error;
206: /* extension -> from local to parallel */
207: VecSet(pcis->vec1_global,0.0);
208: VecSetRandom(pcis->vec1_B,NULL);
209: VecScatterBegin(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
210: VecScatterEnd(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
211: VecDuplicate(pcis->vec1_global,&vec2_global);
212: VecCopy(pcis->vec1_global,vec2_global);
214: VecScatterBegin(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
215: VecScatterEnd(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
216: PCBDDCScalingExtension(pc,pcis->vec1_B,pcis->vec1_global);
217: VecAXPY(pcis->vec1_global,-1.0,vec2_global);
218: VecNorm(pcis->vec1_global,NORM_INFINITY,&error);
219: PetscViewerASCIIPrintf(viewer,"Error scaling extension %1.14e\n",error);
220: if (error>1.e-8 && pcbddc->dbg_flag>1) {
221: VecView(pcis->vec1_global,viewer);
222: }
223: VecDestroy(&vec2_global);
225: /* restriction -> from parallel to local */
226: VecSet(pcis->vec1_global,0.0);
227: VecSetRandom(pcis->vec1_B,NULL);
228: VecScatterBegin(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
229: VecScatterEnd(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
231: PCBDDCScalingRestriction(pc,pcis->vec1_global,pcis->vec1_B);
232: VecScale(pcis->vec1_B,-1.0);
233: VecScatterBegin(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
234: VecScatterEnd(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
235: VecNorm(pcis->vec1_global,NORM_INFINITY,&error);
236: PetscViewerASCIIPrintf(viewer,"Error scaling restriction %1.14e\n",error);
237: if (error>1.e-8 && pcbddc->dbg_flag>1) {
238: VecView(pcis->vec1_global,viewer);
239: }
240: }
241: return(0);
242: }
246: PetscErrorCode PCBDDCScalingDestroy(PC pc)
247: {
248: PC_BDDC* pcbddc=(PC_BDDC*)pc->data;
252: if (pcbddc->deluxe_ctx) {
253: PCBDDCScalingDestroy_Deluxe(pc);
254: }
255: VecDestroy(&pcbddc->work_scaling);
256: /* remove functions */
257: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingRestriction_C",NULL);
258: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingExtension_C",NULL);
259: return(0);
260: }
264: static PetscErrorCode PCBDDCScalingCreate_Deluxe(PC pc)
265: {
266: PC_BDDC* pcbddc=(PC_BDDC*)pc->data;
267: PCBDDCDeluxeScaling deluxe_ctx;
268: PetscErrorCode ierr;
271: PetscNew(&deluxe_ctx);
272: pcbddc->deluxe_ctx = deluxe_ctx;
273: return(0);
274: }
278: static PetscErrorCode PCBDDCScalingDestroy_Deluxe(PC pc)
279: {
280: PC_BDDC* pcbddc=(PC_BDDC*)pc->data;
281: PetscErrorCode ierr;
284: PCBDDCScalingReset_Deluxe_Solvers(pcbddc->deluxe_ctx);
285: PetscFree(pcbddc->deluxe_ctx);
286: return(0);
287: }
291: static PetscErrorCode PCBDDCScalingReset_Deluxe_Solvers(PCBDDCDeluxeScaling deluxe_ctx)
292: {
293: PetscErrorCode ierr;
296: PetscFree(deluxe_ctx->idx_simple_B);
297: deluxe_ctx->n_simple = 0;
298: VecScatterDestroy(&deluxe_ctx->seq_scctx);
299: VecDestroy(&deluxe_ctx->seq_work1);
300: VecDestroy(&deluxe_ctx->seq_work2);
301: MatDestroy(&deluxe_ctx->seq_mat);
302: KSPDestroy(&deluxe_ctx->seq_ksp);
303: return(0);
304: }
308: static PetscErrorCode PCBDDCScalingSetUp_Deluxe(PC pc)
309: {
310: PC_IS *pcis=(PC_IS*)pc->data;
311: PC_BDDC *pcbddc=(PC_BDDC*)pc->data;
312: PCBDDCDeluxeScaling deluxe_ctx=pcbddc->deluxe_ctx;
313: PCBDDCSubSchurs sub_schurs=pcbddc->sub_schurs;
314: PetscErrorCode ierr;
317: /* (TODO: reuse) throw away the solvers */
318: PCBDDCScalingReset_Deluxe_Solvers(deluxe_ctx);
320: /* Compute data structures to solve sequential problems */
321: PCBDDCScalingSetUp_Deluxe_Private(pc);
323: /* diagonal scaling on interface dofs not contained in cc */
324: if (sub_schurs->is_vertices || sub_schurs->is_dir) {
325: PetscInt n_com,n_dir;
326: n_com = 0;
327: if (sub_schurs->is_vertices) {
328: ISGetLocalSize(sub_schurs->is_vertices,&n_com);
329: }
330: n_dir = 0;
331: if (sub_schurs->is_dir) {
332: ISGetLocalSize(sub_schurs->is_dir,&n_dir);
333: }
334: deluxe_ctx->n_simple = n_dir + n_com;
335: PetscMalloc1(deluxe_ctx->n_simple,&deluxe_ctx->idx_simple_B);
336: if (sub_schurs->is_vertices) {
337: PetscInt nmap;
338: const PetscInt *idxs;
340: ISGetIndices(sub_schurs->is_vertices,&idxs);
341: ISGlobalToLocalMappingApply(pcis->BtoNmap,IS_GTOLM_DROP,n_com,idxs,&nmap,deluxe_ctx->idx_simple_B);
342: if (nmap != n_com) {
343: SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error when mapping simply scaled dofs (is_vertices)! %d != %d",nmap,n_com);
344: }
345: ISRestoreIndices(sub_schurs->is_vertices,&idxs);
346: }
347: if (sub_schurs->is_dir) {
348: PetscInt nmap;
349: const PetscInt *idxs;
351: ISGetIndices(sub_schurs->is_dir,&idxs);
352: ISGlobalToLocalMappingApply(pcis->BtoNmap,IS_GTOLM_DROP,n_dir,idxs,&nmap,deluxe_ctx->idx_simple_B+n_com);
353: if (nmap != n_dir) {
354: SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error when mapping simply scaled dofs (sub_schurs->is_dir)! %d != %d",nmap,n_dir);
355: }
356: ISRestoreIndices(sub_schurs->is_dir,&idxs);
357: }
358: PetscSortInt(deluxe_ctx->n_simple,deluxe_ctx->idx_simple_B);
359: } else {
360: deluxe_ctx->n_simple = 0;
361: deluxe_ctx->idx_simple_B = 0;
362: }
363: return(0);
364: }
368: static PetscErrorCode PCBDDCScalingSetUp_Deluxe_Private(PC pc)
369: {
370: PC_BDDC *pcbddc=(PC_BDDC*)pc->data;
371: PCBDDCDeluxeScaling deluxe_ctx=pcbddc->deluxe_ctx;
372: PCBDDCSubSchurs sub_schurs = pcbddc->sub_schurs;
373: PetscErrorCode ierr;
376: if (!sub_schurs->n_subs) {
377: return(0);
378: }
380: /* Create work vectors for sequential part of deluxe */
381: MatCreateVecs(sub_schurs->S_Ej_all,&deluxe_ctx->seq_work1,&deluxe_ctx->seq_work2);
383: /* Compute deluxe sequential scatter */
384: if (sub_schurs->reuse_mumps && !sub_schurs->is_dir) {
385: PCBDDCReuseMumps reuse_mumps = sub_schurs->reuse_mumps;
386: PetscObjectReference((PetscObject)reuse_mumps->correction_scatter_B);
387: deluxe_ctx->seq_scctx = reuse_mumps->correction_scatter_B;
388: } else {
389: VecScatterCreate(pcbddc->work_scaling,sub_schurs->is_Ej_all,deluxe_ctx->seq_work1,NULL,&deluxe_ctx->seq_scctx);
390: }
392: /* Create Mat object for deluxe scaling */
393: PetscObjectReference((PetscObject)sub_schurs->S_Ej_all);
394: deluxe_ctx->seq_mat = sub_schurs->S_Ej_all;
395: if (sub_schurs->sum_S_Ej_all) { /* if this matrix is present, then we need to create the KSP object to invert it */
396: PC pc_temp;
397: MatSolverPackage solver=NULL;
398: char ksp_prefix[256];
399: size_t len;
401: KSPCreate(PETSC_COMM_SELF,&deluxe_ctx->seq_ksp);
402: KSPSetOperators(deluxe_ctx->seq_ksp,sub_schurs->sum_S_Ej_all,sub_schurs->sum_S_Ej_all);
403: KSPSetType(deluxe_ctx->seq_ksp,KSPPREONLY);
404: KSPGetPC(deluxe_ctx->seq_ksp,&pc_temp);
405: PCSetType(pc_temp,PCLU);
406: KSPGetPC(pcbddc->ksp_D,&pc_temp);
407: PCFactorGetMatSolverPackage(pc_temp,(const MatSolverPackage*)&solver);
408: if (solver) {
409: PC new_pc;
410: PCType type;
412: PCGetType(pc_temp,&type);
413: KSPGetPC(deluxe_ctx->seq_ksp,&new_pc);
414: PCSetType(new_pc,type);
415: PCFactorSetMatSolverPackage(new_pc,solver);
416: }
417: PetscStrlen(((PetscObject)(pcbddc->ksp_D))->prefix,&len);
418: len -= 10; /* remove "dirichlet_" */
419: PetscStrncpy(ksp_prefix,((PetscObject)(pcbddc->ksp_D))->prefix,len+1);
420: PetscStrcat(ksp_prefix,"deluxe_");
421: KSPSetOptionsPrefix(deluxe_ctx->seq_ksp,ksp_prefix);
422: KSPSetFromOptions(deluxe_ctx->seq_ksp);
423: KSPSetUp(deluxe_ctx->seq_ksp);
424: }
425: return(0);
426: }