Actual source code: bddcscalingbasic.c
petsc-3.7.3 2016-08-01
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) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Local vector cannot be pcbddc->work_scaling!\n");
83: PetscUseMethod(pc,"PCBDDCScalingExtension_C",(PC,Vec,Vec),(pc,local_interface_vector,global_vector));
84: return(0);
85: }
89: static PetscErrorCode PCBDDCScalingRestriction_Basic(PC pc, Vec global_vector, Vec local_interface_vector)
90: {
92: PC_IS *pcis = (PC_IS*)pc->data;
95: VecScatterBegin(pcis->global_to_B,global_vector,local_interface_vector,INSERT_VALUES,SCATTER_FORWARD);
96: VecScatterEnd(pcis->global_to_B,global_vector,local_interface_vector,INSERT_VALUES,SCATTER_FORWARD);
97: /* Apply partition of unity */
98: VecPointwiseMult(local_interface_vector,pcis->D,local_interface_vector);
99: return(0);
100: }
104: static PetscErrorCode PCBDDCScalingRestriction_Deluxe(PC pc, Vec x, Vec y)
105: {
106: PC_IS* pcis=(PC_IS*)pc->data;
107: PC_BDDC* pcbddc=(PC_BDDC*)pc->data;
108: PCBDDCDeluxeScaling deluxe_ctx = pcbddc->deluxe_ctx;
109: PetscErrorCode ierr;
112: /* get local boundary part of global vector */
113: VecScatterBegin(pcis->global_to_B,x,y,INSERT_VALUES,SCATTER_FORWARD);
114: VecScatterEnd(pcis->global_to_B,x,y,INSERT_VALUES,SCATTER_FORWARD);
115: if (deluxe_ctx->n_simple) { /* scale deluxe vertices using diagonal scaling */
116: PetscInt i;
117: PetscScalar *array_y;
118: const PetscScalar *array_D;
119: VecGetArray(y,&array_y);
120: VecGetArrayRead(pcis->D,&array_D);
121: for (i=0;i<deluxe_ctx->n_simple;i++) {
122: array_y[deluxe_ctx->idx_simple_B[i]] *= array_D[deluxe_ctx->idx_simple_B[i]];
123: }
124: VecRestoreArrayRead(pcis->D,&array_D);
125: VecRestoreArray(y,&array_y);
126: }
127: /* sequential part : all problems and Schur applications collapsed into a single matrix vector multiplication and ksp solution */
128: if (deluxe_ctx->seq_mat) {
129: VecScatterBegin(deluxe_ctx->seq_scctx,y,deluxe_ctx->seq_work1,INSERT_VALUES,SCATTER_FORWARD);
130: VecScatterEnd(deluxe_ctx->seq_scctx,y,deluxe_ctx->seq_work1,INSERT_VALUES,SCATTER_FORWARD);
131: if (deluxe_ctx->seq_ksp) {
132: KSPSolve(deluxe_ctx->seq_ksp,deluxe_ctx->seq_work1,deluxe_ctx->seq_work1);
133: }
134: MatMult(deluxe_ctx->seq_mat,deluxe_ctx->seq_work1,deluxe_ctx->seq_work2);
135: VecScatterBegin(deluxe_ctx->seq_scctx,deluxe_ctx->seq_work2,y,INSERT_VALUES,SCATTER_REVERSE);
136: VecScatterEnd(deluxe_ctx->seq_scctx,deluxe_ctx->seq_work2,y,INSERT_VALUES,SCATTER_REVERSE);
137: }
138: return(0);
139: }
143: PetscErrorCode PCBDDCScalingRestriction(PC pc, Vec global_vector, Vec local_interface_vector)
144: {
145: PC_BDDC *pcbddc=(PC_BDDC*)pc->data;
152: if (local_interface_vector == pcbddc->work_scaling) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Local vector cannot be pcbddc->work_scaling!\n");
153: PetscUseMethod(pc,"PCBDDCScalingRestriction_C",(PC,Vec,Vec),(pc,global_vector,local_interface_vector));
154: return(0);
155: }
159: PetscErrorCode PCBDDCScalingSetUp(PC pc)
160: {
161: PC_IS* pcis=(PC_IS*)pc->data;
162: PC_BDDC* pcbddc=(PC_BDDC*)pc->data;
167: /* create work vector for the operator */
168: VecDestroy(&pcbddc->work_scaling);
169: VecDuplicate(pcis->vec1_B,&pcbddc->work_scaling);
170: /* always rebuild pcis->D */
171: if (pcis->use_stiffness_scaling) {
172: MatGetDiagonal(pcbddc->local_mat,pcis->vec1_N);
173: VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);
174: VecScatterEnd(pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);
175: }
176: VecCopy(pcis->D,pcis->vec1_B);
177: VecSet(pcis->vec1_global,0.0);
178: VecScatterBegin(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
179: VecScatterEnd(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
180: VecScatterBegin(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
181: VecScatterEnd(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
182: VecPointwiseDivide(pcis->D,pcis->D,pcis->vec1_B);
183: /* now setup */
184: if (pcbddc->use_deluxe_scaling) {
185: if (!pcbddc->deluxe_ctx) {
186: PCBDDCScalingCreate_Deluxe(pc);
187: }
188: PCBDDCScalingSetUp_Deluxe(pc);
189: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingRestriction_C",PCBDDCScalingRestriction_Deluxe);
190: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingExtension_C",PCBDDCScalingExtension_Deluxe);
191: } else {
192: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingRestriction_C",PCBDDCScalingRestriction_Basic);
193: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingExtension_C",PCBDDCScalingExtension_Basic);
194: }
196: /* test */
197: if (pcbddc->dbg_flag) {
198: Vec vec2_global;
199: PetscViewer viewer=pcbddc->dbg_viewer;
200: PetscReal error;
202: /* extension -> from local to parallel */
203: VecSet(pcis->vec1_global,0.0);
204: VecSetRandom(pcis->vec1_B,NULL);
205: VecScatterBegin(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
206: VecScatterEnd(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
207: VecDuplicate(pcis->vec1_global,&vec2_global);
208: VecCopy(pcis->vec1_global,vec2_global);
210: VecScatterBegin(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
211: VecScatterEnd(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
212: PCBDDCScalingExtension(pc,pcis->vec1_B,pcis->vec1_global);
213: VecAXPY(pcis->vec1_global,-1.0,vec2_global);
214: VecNorm(pcis->vec1_global,NORM_INFINITY,&error);
215: PetscViewerASCIIPrintf(viewer,"Error scaling extension %1.14e\n",error);
216: if (error>1.e-8 && pcbddc->dbg_flag>1) {
217: VecView(pcis->vec1_global,viewer);
218: }
219: VecDestroy(&vec2_global);
221: /* restriction -> from parallel to local */
222: VecSet(pcis->vec1_global,0.0);
223: VecSetRandom(pcis->vec1_B,NULL);
224: VecScatterBegin(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
225: VecScatterEnd(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
227: PCBDDCScalingRestriction(pc,pcis->vec1_global,pcis->vec1_B);
228: VecScale(pcis->vec1_B,-1.0);
229: VecScatterBegin(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
230: VecScatterEnd(pcis->global_to_B,pcis->vec1_B,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
231: VecNorm(pcis->vec1_global,NORM_INFINITY,&error);
232: PetscViewerASCIIPrintf(viewer,"Error scaling restriction %1.14e\n",error);
233: if (error>1.e-8 && pcbddc->dbg_flag>1) {
234: VecView(pcis->vec1_global,viewer);
235: }
236: }
237: return(0);
238: }
242: PetscErrorCode PCBDDCScalingDestroy(PC pc)
243: {
244: PC_BDDC* pcbddc=(PC_BDDC*)pc->data;
248: if (pcbddc->deluxe_ctx) {
249: PCBDDCScalingDestroy_Deluxe(pc);
250: }
251: VecDestroy(&pcbddc->work_scaling);
252: /* remove functions */
253: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingRestriction_C",NULL);
254: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCScalingExtension_C",NULL);
255: return(0);
256: }
260: static PetscErrorCode PCBDDCScalingCreate_Deluxe(PC pc)
261: {
262: PC_BDDC* pcbddc=(PC_BDDC*)pc->data;
263: PCBDDCDeluxeScaling deluxe_ctx;
264: PetscErrorCode ierr;
267: PetscNew(&deluxe_ctx);
268: pcbddc->deluxe_ctx = deluxe_ctx;
269: return(0);
270: }
274: static PetscErrorCode PCBDDCScalingDestroy_Deluxe(PC pc)
275: {
276: PC_BDDC* pcbddc=(PC_BDDC*)pc->data;
277: PetscErrorCode ierr;
280: PCBDDCScalingReset_Deluxe_Solvers(pcbddc->deluxe_ctx);
281: PetscFree(pcbddc->deluxe_ctx);
282: return(0);
283: }
287: static PetscErrorCode PCBDDCScalingReset_Deluxe_Solvers(PCBDDCDeluxeScaling deluxe_ctx)
288: {
289: PetscErrorCode ierr;
292: PetscFree(deluxe_ctx->idx_simple_B);
293: deluxe_ctx->n_simple = 0;
294: VecScatterDestroy(&deluxe_ctx->seq_scctx);
295: VecDestroy(&deluxe_ctx->seq_work1);
296: VecDestroy(&deluxe_ctx->seq_work2);
297: MatDestroy(&deluxe_ctx->seq_mat);
298: KSPDestroy(&deluxe_ctx->seq_ksp);
299: return(0);
300: }
304: static PetscErrorCode PCBDDCScalingSetUp_Deluxe(PC pc)
305: {
306: PC_IS *pcis=(PC_IS*)pc->data;
307: PC_BDDC *pcbddc=(PC_BDDC*)pc->data;
308: PCBDDCDeluxeScaling deluxe_ctx=pcbddc->deluxe_ctx;
309: PCBDDCSubSchurs sub_schurs=pcbddc->sub_schurs;
310: PetscErrorCode ierr;
313: /* (TODO: reuse) throw away the solvers */
314: PCBDDCScalingReset_Deluxe_Solvers(deluxe_ctx);
316: /* Compute data structures to solve sequential problems */
317: PCBDDCScalingSetUp_Deluxe_Private(pc);
319: /* diagonal scaling on interface dofs not contained in cc */
320: if (sub_schurs->is_vertices || sub_schurs->is_dir) {
321: PetscInt n_com,n_dir;
322: n_com = 0;
323: if (sub_schurs->is_vertices) {
324: ISGetLocalSize(sub_schurs->is_vertices,&n_com);
325: }
326: n_dir = 0;
327: if (sub_schurs->is_dir) {
328: ISGetLocalSize(sub_schurs->is_dir,&n_dir);
329: }
330: deluxe_ctx->n_simple = n_dir + n_com;
331: PetscMalloc1(deluxe_ctx->n_simple,&deluxe_ctx->idx_simple_B);
332: if (sub_schurs->is_vertices) {
333: PetscInt nmap;
334: const PetscInt *idxs;
336: ISGetIndices(sub_schurs->is_vertices,&idxs);
337: ISGlobalToLocalMappingApply(pcis->BtoNmap,IS_GTOLM_DROP,n_com,idxs,&nmap,deluxe_ctx->idx_simple_B);
338: if (nmap != n_com) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error when mapping simply scaled dofs (is_vertices)! %D != %D",nmap,n_com);
339: ISRestoreIndices(sub_schurs->is_vertices,&idxs);
340: }
341: if (sub_schurs->is_dir) {
342: PetscInt nmap;
343: const PetscInt *idxs;
345: ISGetIndices(sub_schurs->is_dir,&idxs);
346: ISGlobalToLocalMappingApply(pcis->BtoNmap,IS_GTOLM_DROP,n_dir,idxs,&nmap,deluxe_ctx->idx_simple_B+n_com);
347: 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);
348: ISRestoreIndices(sub_schurs->is_dir,&idxs);
349: }
350: PetscSortInt(deluxe_ctx->n_simple,deluxe_ctx->idx_simple_B);
351: } else {
352: deluxe_ctx->n_simple = 0;
353: deluxe_ctx->idx_simple_B = 0;
354: }
355: return(0);
356: }
360: static PetscErrorCode PCBDDCScalingSetUp_Deluxe_Private(PC pc)
361: {
362: PC_BDDC *pcbddc=(PC_BDDC*)pc->data;
363: PCBDDCDeluxeScaling deluxe_ctx=pcbddc->deluxe_ctx;
364: PCBDDCSubSchurs sub_schurs = pcbddc->sub_schurs;
365: PetscErrorCode ierr;
368: if (!sub_schurs->n_subs) {
369: return(0);
370: }
372: /* Create work vectors for sequential part of deluxe */
373: MatCreateVecs(sub_schurs->S_Ej_all,&deluxe_ctx->seq_work1,&deluxe_ctx->seq_work2);
375: /* Compute deluxe sequential scatter */
376: if (sub_schurs->reuse_mumps && !sub_schurs->is_dir) {
377: PCBDDCReuseMumps reuse_mumps = sub_schurs->reuse_mumps;
378: PetscObjectReference((PetscObject)reuse_mumps->correction_scatter_B);
379: deluxe_ctx->seq_scctx = reuse_mumps->correction_scatter_B;
380: } else {
381: VecScatterCreate(pcbddc->work_scaling,sub_schurs->is_Ej_all,deluxe_ctx->seq_work1,NULL,&deluxe_ctx->seq_scctx);
382: }
384: /* Create Mat object for deluxe scaling */
385: PetscObjectReference((PetscObject)sub_schurs->S_Ej_all);
386: deluxe_ctx->seq_mat = sub_schurs->S_Ej_all;
387: if (sub_schurs->sum_S_Ej_all) { /* if this matrix is present, then we need to create the KSP object to invert it */
388: PC pc_temp;
389: MatSolverPackage solver=NULL;
390: char ksp_prefix[256];
391: size_t len;
393: KSPCreate(PETSC_COMM_SELF,&deluxe_ctx->seq_ksp);
394: KSPSetOperators(deluxe_ctx->seq_ksp,sub_schurs->sum_S_Ej_all,sub_schurs->sum_S_Ej_all);
395: KSPSetType(deluxe_ctx->seq_ksp,KSPPREONLY);
396: KSPGetPC(deluxe_ctx->seq_ksp,&pc_temp);
397: PCSetType(pc_temp,PCLU);
398: KSPGetPC(pcbddc->ksp_D,&pc_temp);
399: PCFactorGetMatSolverPackage(pc_temp,(const MatSolverPackage*)&solver);
400: if (solver) {
401: PC new_pc;
402: PCType type;
404: PCGetType(pc_temp,&type);
405: KSPGetPC(deluxe_ctx->seq_ksp,&new_pc);
406: PCSetType(new_pc,type);
407: PCFactorSetMatSolverPackage(new_pc,solver);
408: }
409: PetscStrlen(((PetscObject)(pcbddc->ksp_D))->prefix,&len);
410: len -= 10; /* remove "dirichlet_" */
411: PetscStrncpy(ksp_prefix,((PetscObject)(pcbddc->ksp_D))->prefix,len+1);
412: PetscStrcat(ksp_prefix,"deluxe_");
413: KSPSetOptionsPrefix(deluxe_ctx->seq_ksp,ksp_prefix);
414: KSPSetFromOptions(deluxe_ctx->seq_ksp);
415: KSPSetUp(deluxe_ctx->seq_ksp);
416: }
417: return(0);
418: }