Actual source code: redundant.c
petsc-3.6.4 2016-04-12
2: /*
3: This file defines a "solve the problem redundantly on each subgroup of processor" preconditioner.
4: */
5: #include <petsc/private/pcimpl.h>
6: #include <petscksp.h> /*I "petscksp.h" I*/
8: typedef struct {
9: KSP ksp;
10: PC pc; /* actual preconditioner used on each processor */
11: Vec xsub,ysub; /* vectors of a subcommunicator to hold parallel vectors of PetscObjectComm((PetscObject)pc) */
12: Vec xdup,ydup; /* parallel vector that congregates xsub or ysub facilitating vector scattering */
13: Mat pmats; /* matrix and optional preconditioner matrix belong to a subcommunicator */
14: VecScatter scatterin,scatterout; /* scatter used to move all values to each processor group (subcommunicator) */
15: PetscBool useparallelmat;
16: PetscSubcomm psubcomm;
17: PetscInt nsubcomm; /* num of data structure PetscSubcomm */
18: } PC_Redundant;
22: static PetscErrorCode PCView_Redundant(PC pc,PetscViewer viewer)
23: {
24: PC_Redundant *red = (PC_Redundant*)pc->data;
26: PetscBool iascii,isstring;
27: PetscViewer subviewer;
30: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
31: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
32: if (iascii) {
33: if (!red->psubcomm) {
34: PetscViewerASCIIPrintf(viewer," Redundant preconditioner: Not yet setup\n");
35: } else {
36: PetscViewerASCIIPrintf(viewer," Redundant preconditioner: First (color=0) of %D PCs follows\n",red->nsubcomm);
37: PetscViewerGetSubcomm(viewer,((PetscObject)red->pc)->comm,&subviewer);
38: if (!red->psubcomm->color) { /* only view first redundant pc */
39: PetscViewerASCIIPushTab(viewer);
40: KSPView(red->ksp,subviewer);
41: PetscViewerASCIIPopTab(viewer);
42: }
43: PetscViewerRestoreSubcomm(viewer,((PetscObject)red->pc)->comm,&subviewer);
44: }
45: } else if (isstring) {
46: PetscViewerStringSPrintf(viewer," Redundant solver preconditioner");
47: }
48: return(0);
49: }
53: static PetscErrorCode PCSetUp_Redundant(PC pc)
54: {
55: PC_Redundant *red = (PC_Redundant*)pc->data;
57: PetscInt mstart,mend,mlocal,M;
58: PetscMPIInt size;
59: MPI_Comm comm,subcomm;
60: Vec x;
61: const char *prefix;
64: PetscObjectGetComm((PetscObject)pc,&comm);
65:
66: /* if pmatrix set by user is sequential then we do not need to gather the parallel matrix */
67: MPI_Comm_size(comm,&size);
68: if (size == 1) red->useparallelmat = PETSC_FALSE;
70: if (!pc->setupcalled) {
71: PetscInt mloc_sub;
72: if (!red->psubcomm) {
73: PetscSubcommCreate(comm,&red->psubcomm);
74: PetscSubcommSetNumber(red->psubcomm,red->nsubcomm);
75: PetscSubcommSetType(red->psubcomm,PETSC_SUBCOMM_CONTIGUOUS);
76: /* enable runtime switch of psubcomm type, e.g., '-psubcomm_type interlaced */
77: PetscSubcommSetFromOptions(red->psubcomm);
78: PetscLogObjectMemory((PetscObject)pc,sizeof(PetscSubcomm));
80: /* create a new PC that processors in each subcomm have copy of */
81: subcomm = PetscSubcommChild(red->psubcomm);
83: KSPCreate(subcomm,&red->ksp);
84: KSPSetErrorIfNotConverged(red->ksp,pc->erroriffailure);
85: PetscObjectIncrementTabLevel((PetscObject)red->ksp,(PetscObject)pc,1);
86: PetscLogObjectParent((PetscObject)pc,(PetscObject)red->ksp);
87: KSPSetType(red->ksp,KSPPREONLY);
88: KSPGetPC(red->ksp,&red->pc);
89: PCSetType(red->pc,PCLU);
91: PCGetOptionsPrefix(pc,&prefix);
92: KSPSetOptionsPrefix(red->ksp,prefix);
93: KSPAppendOptionsPrefix(red->ksp,"redundant_");
94: } else {
95: subcomm = PetscSubcommChild(red->psubcomm);
96: }
98: if (red->useparallelmat) {
99: /* grab the parallel matrix and put it into processors of a subcomminicator */
100: MatCreateRedundantMatrix(pc->pmat,red->psubcomm->n,subcomm,MAT_INITIAL_MATRIX,&red->pmats);
101: KSPSetOperators(red->ksp,red->pmats,red->pmats);
102:
103: /* get working vectors xsub and ysub */
104: MatCreateVecs(red->pmats,&red->xsub,&red->ysub);
106: /* create working vectors xdup and ydup.
107: xdup concatenates all xsub's contigously to form a mpi vector over dupcomm (see PetscSubcommCreate_interlaced())
108: ydup concatenates all ysub and has empty local arrays because ysub's arrays will be place into it.
109: Note: we use communicator dupcomm, not PetscObjectComm((PetscObject)pc)! */
110: MatGetLocalSize(red->pmats,&mloc_sub,NULL);
111: VecCreateMPI(PetscSubcommContiguousParent(red->psubcomm),mloc_sub,PETSC_DECIDE,&red->xdup);
112: VecCreateMPIWithArray(PetscSubcommContiguousParent(red->psubcomm),1,mloc_sub,PETSC_DECIDE,NULL,&red->ydup);
114: /* create vecscatters */
115: if (!red->scatterin) { /* efficiency of scatterin is independent from psubcomm_type! */
116: IS is1,is2;
117: PetscInt *idx1,*idx2,i,j,k;
119: MatCreateVecs(pc->pmat,&x,0);
120: VecGetSize(x,&M);
121: VecGetOwnershipRange(x,&mstart,&mend);
122: mlocal = mend - mstart;
123: PetscMalloc2(red->psubcomm->n*mlocal,&idx1,red->psubcomm->n*mlocal,&idx2);
124: j = 0;
125: for (k=0; k<red->psubcomm->n; k++) {
126: for (i=mstart; i<mend; i++) {
127: idx1[j] = i;
128: idx2[j++] = i + M*k;
129: }
130: }
131: ISCreateGeneral(comm,red->psubcomm->n*mlocal,idx1,PETSC_COPY_VALUES,&is1);
132: ISCreateGeneral(comm,red->psubcomm->n*mlocal,idx2,PETSC_COPY_VALUES,&is2);
133: VecScatterCreate(x,is1,red->xdup,is2,&red->scatterin);
134: ISDestroy(&is1);
135: ISDestroy(&is2);
137: /* Impl below is good for PETSC_SUBCOMM_INTERLACED (no inter-process communication) and PETSC_SUBCOMM_CONTIGUOUS (communication within subcomm) */
138: ISCreateStride(comm,mlocal,mstart+ red->psubcomm->color*M,1,&is1);
139: ISCreateStride(comm,mlocal,mstart,1,&is2);
140: VecScatterCreate(red->xdup,is1,x,is2,&red->scatterout);
141: ISDestroy(&is1);
142: ISDestroy(&is2);
143: PetscFree2(idx1,idx2);
144: VecDestroy(&x);
145: }
146: } else { /* !red->useparallelmat */
147: KSPSetOperators(red->ksp,pc->mat,pc->pmat);
148: }
149: } else { /* pc->setupcalled */
150: if (red->useparallelmat) {
151: MatReuse reuse;
152: /* grab the parallel matrix and put it into processors of a subcomminicator */
153: /*--------------------------------------------------------------------------*/
154: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
155: /* destroy old matrices */
156: MatDestroy(&red->pmats);
157: reuse = MAT_INITIAL_MATRIX;
158: } else {
159: reuse = MAT_REUSE_MATRIX;
160: }
161: MatCreateRedundantMatrix(pc->pmat,red->psubcomm->n,PetscSubcommChild(red->psubcomm),reuse,&red->pmats);
162: KSPSetOperators(red->ksp,red->pmats,red->pmats);
163: } else { /* !red->useparallelmat */
164: KSPSetOperators(red->ksp,pc->mat,pc->pmat);
165: }
166: }
168: if (pc->setfromoptionscalled) {
169: KSPSetFromOptions(red->ksp);
170: }
171: KSPSetUp(red->ksp);
172: return(0);
173: }
177: static PetscErrorCode PCApply_Redundant(PC pc,Vec x,Vec y)
178: {
179: PC_Redundant *red = (PC_Redundant*)pc->data;
181: PetscScalar *array;
184: if (!red->useparallelmat) {
185: KSPSolve(red->ksp,x,y);
186: return(0);
187: }
189: /* scatter x to xdup */
190: VecScatterBegin(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD);
191: VecScatterEnd(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD);
193: /* place xdup's local array into xsub */
194: VecGetArray(red->xdup,&array);
195: VecPlaceArray(red->xsub,(const PetscScalar*)array);
197: /* apply preconditioner on each processor */
198: KSPSolve(red->ksp,red->xsub,red->ysub);
199: VecResetArray(red->xsub);
200: VecRestoreArray(red->xdup,&array);
202: /* place ysub's local array into ydup */
203: VecGetArray(red->ysub,&array);
204: VecPlaceArray(red->ydup,(const PetscScalar*)array);
206: /* scatter ydup to y */
207: VecScatterBegin(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD);
208: VecScatterEnd(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD);
209: VecResetArray(red->ydup);
210: VecRestoreArray(red->ysub,&array);
211: return(0);
212: }
216: static PetscErrorCode PCApplyTranspose_Redundant(PC pc,Vec x,Vec y)
217: {
218: PC_Redundant *red = (PC_Redundant*)pc->data;
220: PetscScalar *array;
223: if (!red->useparallelmat) {
224: KSPSolveTranspose(red->ksp,x,y);
225: return(0);
226: }
228: /* scatter x to xdup */
229: VecScatterBegin(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD);
230: VecScatterEnd(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD);
232: /* place xdup's local array into xsub */
233: VecGetArray(red->xdup,&array);
234: VecPlaceArray(red->xsub,(const PetscScalar*)array);
236: /* apply preconditioner on each processor */
237: KSPSolveTranspose(red->ksp,red->xsub,red->ysub);
238: VecResetArray(red->xsub);
239: VecRestoreArray(red->xdup,&array);
241: /* place ysub's local array into ydup */
242: VecGetArray(red->ysub,&array);
243: VecPlaceArray(red->ydup,(const PetscScalar*)array);
245: /* scatter ydup to y */
246: VecScatterBegin(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD);
247: VecScatterEnd(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD);
248: VecResetArray(red->ydup);
249: VecRestoreArray(red->ysub,&array);
250: return(0);
251: }
255: static PetscErrorCode PCReset_Redundant(PC pc)
256: {
257: PC_Redundant *red = (PC_Redundant*)pc->data;
261: if (red->useparallelmat) {
262: VecScatterDestroy(&red->scatterin);
263: VecScatterDestroy(&red->scatterout);
264: VecDestroy(&red->ysub);
265: VecDestroy(&red->xsub);
266: VecDestroy(&red->xdup);
267: VecDestroy(&red->ydup);
268: }
269: MatDestroy(&red->pmats);
270: KSPReset(red->ksp);
271: return(0);
272: }
276: static PetscErrorCode PCDestroy_Redundant(PC pc)
277: {
278: PC_Redundant *red = (PC_Redundant*)pc->data;
282: PCReset_Redundant(pc);
283: KSPDestroy(&red->ksp);
284: PetscSubcommDestroy(&red->psubcomm);
285: PetscFree(pc->data);
286: return(0);
287: }
291: static PetscErrorCode PCSetFromOptions_Redundant(PetscOptions *PetscOptionsObject,PC pc)
292: {
294: PC_Redundant *red = (PC_Redundant*)pc->data;
297: PetscOptionsHead(PetscOptionsObject,"Redundant options");
298: PetscOptionsInt("-pc_redundant_number","Number of redundant pc","PCRedundantSetNumber",red->nsubcomm,&red->nsubcomm,0);
299: PetscOptionsTail();
300: return(0);
301: }
305: static PetscErrorCode PCRedundantSetNumber_Redundant(PC pc,PetscInt nreds)
306: {
307: PC_Redundant *red = (PC_Redundant*)pc->data;
310: red->nsubcomm = nreds;
311: return(0);
312: }
316: /*@
317: PCRedundantSetNumber - Sets the number of redundant preconditioner contexts.
319: Logically Collective on PC
321: Input Parameters:
322: + pc - the preconditioner context
323: - nredundant - number of redundant preconditioner contexts; for example if you are using 64 MPI processes and
324: use an nredundant of 4 there will be 4 parallel solves each on 16 = 64/4 processes.
326: Level: advanced
328: .keywords: PC, redundant solve
329: @*/
330: PetscErrorCode PCRedundantSetNumber(PC pc,PetscInt nredundant)
331: {
336: if (nredundant <= 0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "num of redundant pc %D must be positive",nredundant);
337: PetscTryMethod(pc,"PCRedundantSetNumber_C",(PC,PetscInt),(pc,nredundant));
338: return(0);
339: }
343: static PetscErrorCode PCRedundantSetScatter_Redundant(PC pc,VecScatter in,VecScatter out)
344: {
345: PC_Redundant *red = (PC_Redundant*)pc->data;
349: PetscObjectReference((PetscObject)in);
350: VecScatterDestroy(&red->scatterin);
352: red->scatterin = in;
354: PetscObjectReference((PetscObject)out);
355: VecScatterDestroy(&red->scatterout);
356: red->scatterout = out;
357: return(0);
358: }
362: /*@
363: PCRedundantSetScatter - Sets the scatter used to copy values into the
364: redundant local solve and the scatter to move them back into the global
365: vector.
367: Logically Collective on PC
369: Input Parameters:
370: + pc - the preconditioner context
371: . in - the scatter to move the values in
372: - out - the scatter to move them out
374: Level: advanced
376: .keywords: PC, redundant solve
377: @*/
378: PetscErrorCode PCRedundantSetScatter(PC pc,VecScatter in,VecScatter out)
379: {
386: PetscTryMethod(pc,"PCRedundantSetScatter_C",(PC,VecScatter,VecScatter),(pc,in,out));
387: return(0);
388: }
392: static PetscErrorCode PCRedundantGetKSP_Redundant(PC pc,KSP *innerksp)
393: {
395: PC_Redundant *red = (PC_Redundant*)pc->data;
396: MPI_Comm comm,subcomm;
397: const char *prefix;
400: if (!red->psubcomm) {
401: PetscObjectGetComm((PetscObject)pc,&comm);
402: PetscSubcommCreate(comm,&red->psubcomm);
403: PetscSubcommSetNumber(red->psubcomm,red->nsubcomm);
404: PetscSubcommSetType(red->psubcomm,PETSC_SUBCOMM_INTERLACED);
405: PetscLogObjectMemory((PetscObject)pc,sizeof(PetscSubcomm));
407: /* create a new PC that processors in each subcomm have copy of */
408: subcomm = PetscSubcommChild(red->psubcomm);
410: KSPCreate(subcomm,&red->ksp);
411: KSPSetErrorIfNotConverged(red->ksp,pc->erroriffailure);
412: PetscObjectIncrementTabLevel((PetscObject)red->ksp,(PetscObject)pc,1);
413: PetscLogObjectParent((PetscObject)pc,(PetscObject)red->ksp);
414: KSPSetType(red->ksp,KSPPREONLY);
415: KSPGetPC(red->ksp,&red->pc);
416: PCSetType(red->pc,PCLU);
418: PCGetOptionsPrefix(pc,&prefix);
419: KSPSetOptionsPrefix(red->ksp,prefix);
420: KSPAppendOptionsPrefix(red->ksp,"redundant_");
421: }
422: *innerksp = red->ksp;
423: return(0);
424: }
428: /*@
429: PCRedundantGetKSP - Gets the less parallel KSP created by the redundant PC.
431: Not Collective
433: Input Parameter:
434: . pc - the preconditioner context
436: Output Parameter:
437: . innerksp - the KSP on the smaller set of processes
439: Level: advanced
441: .keywords: PC, redundant solve
442: @*/
443: PetscErrorCode PCRedundantGetKSP(PC pc,KSP *innerksp)
444: {
450: PetscTryMethod(pc,"PCRedundantGetKSP_C",(PC,KSP*),(pc,innerksp));
451: return(0);
452: }
456: static PetscErrorCode PCRedundantGetOperators_Redundant(PC pc,Mat *mat,Mat *pmat)
457: {
458: PC_Redundant *red = (PC_Redundant*)pc->data;
461: if (mat) *mat = red->pmats;
462: if (pmat) *pmat = red->pmats;
463: return(0);
464: }
468: /*@
469: PCRedundantGetOperators - gets the sequential matrix and preconditioner matrix
471: Not Collective
473: Input Parameter:
474: . pc - the preconditioner context
476: Output Parameters:
477: + mat - the matrix
478: - pmat - the (possibly different) preconditioner matrix
480: Level: advanced
482: .keywords: PC, redundant solve
483: @*/
484: PetscErrorCode PCRedundantGetOperators(PC pc,Mat *mat,Mat *pmat)
485: {
492: PetscTryMethod(pc,"PCRedundantGetOperators_C",(PC,Mat*,Mat*),(pc,mat,pmat));
493: return(0);
494: }
496: /* -------------------------------------------------------------------------------------*/
497: /*MC
498: PCREDUNDANT - Runs a KSP solver with preconditioner for the entire problem on subgroups of processors
500: Options for the redundant preconditioners can be set with -redundant_pc_xxx for the redundant KSP with -redundant_ksp_xxx
502: Options Database:
503: . -pc_redundant_number <n> - number of redundant solves, for example if you are using 64 MPI processes and
504: use an n of 4 there will be 4 parallel solves each on 16 = 64/4 processes.
506: Level: intermediate
508: Notes: The default KSP is preonly and the default PC is LU.
510: Developer Notes: Note that PCSetInitialGuessNonzero() is not used by this class but likely should be.
512: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCRedundantSetScatter(),
513: PCRedundantGetKSP(), PCRedundantGetOperators(), PCRedundantSetNumber()
514: M*/
518: PETSC_EXTERN PetscErrorCode PCCreate_Redundant(PC pc)
519: {
521: PC_Redundant *red;
522: PetscMPIInt size;
525: PetscNewLog(pc,&red);
526: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
528: red->nsubcomm = size;
529: red->useparallelmat = PETSC_TRUE;
530: pc->data = (void*)red;
532: pc->ops->apply = PCApply_Redundant;
533: pc->ops->applytranspose = PCApplyTranspose_Redundant;
534: pc->ops->setup = PCSetUp_Redundant;
535: pc->ops->destroy = PCDestroy_Redundant;
536: pc->ops->reset = PCReset_Redundant;
537: pc->ops->setfromoptions = PCSetFromOptions_Redundant;
538: pc->ops->view = PCView_Redundant;
540: PetscObjectComposeFunction((PetscObject)pc,"PCRedundantSetScatter_C",PCRedundantSetScatter_Redundant);
541: PetscObjectComposeFunction((PetscObject)pc,"PCRedundantSetNumber_C",PCRedundantSetNumber_Redundant);
542: PetscObjectComposeFunction((PetscObject)pc,"PCRedundantGetKSP_C",PCRedundantGetKSP_Redundant);
543: PetscObjectComposeFunction((PetscObject)pc,"PCRedundantGetOperators_C",PCRedundantGetOperators_Redundant);
544: return(0);
545: }