Actual source code: redundant.c
petsc-3.4.5 2014-06-29
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,mlocal_sub,rstart_sub,rend_sub,mloc_sub;
58: PetscMPIInt size;
59: MatReuse reuse = MAT_INITIAL_MATRIX;
60: MatStructure str = DIFFERENT_NONZERO_PATTERN;
61: MPI_Comm comm,subcomm;
62: Vec vec;
63: PetscMPIInt subsize,subrank;
64: const char *prefix;
65: const PetscInt *range;
68: PetscObjectGetComm((PetscObject)pc,&comm);
69: MatGetVecs(pc->pmat,&vec,0);
70: VecGetSize(vec,&m);
72: if (!pc->setupcalled) {
73: if (!red->psubcomm) {
74: PetscSubcommCreate(comm,&red->psubcomm);
75: PetscSubcommSetNumber(red->psubcomm,red->nsubcomm);
76: PetscSubcommSetType(red->psubcomm,PETSC_SUBCOMM_INTERLACED);
77: PetscLogObjectMemory(pc,sizeof(PetscSubcomm));
79: /* create a new PC that processors in each subcomm have copy of */
80: subcomm = red->psubcomm->comm;
82: KSPCreate(subcomm,&red->ksp);
83: PetscObjectIncrementTabLevel((PetscObject)red->ksp,(PetscObject)pc,1);
84: PetscLogObjectParent(pc,red->ksp);
85: KSPSetType(red->ksp,KSPPREONLY);
86: KSPGetPC(red->ksp,&red->pc);
87: PCSetType(red->pc,PCLU);
89: PCGetOptionsPrefix(pc,&prefix);
90: KSPSetOptionsPrefix(red->ksp,prefix);
91: KSPAppendOptionsPrefix(red->ksp,"redundant_");
92: } else subcomm = red->psubcomm->comm;
94: /* create working vectors xsub/ysub and xdup/ydup */
95: VecGetLocalSize(vec,&mlocal);
96: VecGetOwnershipRange(vec,&mstart,&mend);
98: /* get local size of xsub/ysub */
99: MPI_Comm_size(subcomm,&subsize);
100: MPI_Comm_rank(subcomm,&subrank);
101: MatGetOwnershipRanges(pc->pmat,&range);
102: rstart_sub = range[red->psubcomm->n*subrank]; /* rstart in xsub/ysub */
103: if (subrank+1 < subsize) rend_sub = range[red->psubcomm->n*(subrank+1)];
104: else rend_sub = m;
106: mloc_sub = rend_sub - rstart_sub;
107: VecCreateMPI(subcomm,mloc_sub,PETSC_DECIDE,&red->ysub);
108: /* create xsub with empty local arrays, because xdup's arrays will be placed into it */
109: VecCreateMPIWithArray(subcomm,1,mloc_sub,PETSC_DECIDE,NULL,&red->xsub);
111: /* create xdup and ydup. ydup has empty local arrays because ysub's arrays will be place into it.
112: Note: we use communicator dupcomm, not PetscObjectComm((PetscObject)pc)! */
113: VecCreateMPI(red->psubcomm->dupparent,mloc_sub,PETSC_DECIDE,&red->xdup);
114: VecCreateMPIWithArray(red->psubcomm->dupparent,1,mloc_sub,PETSC_DECIDE,NULL,&red->ydup);
116: /* create vec scatters */
117: if (!red->scatterin) {
118: IS is1,is2;
119: PetscInt *idx1,*idx2,i,j,k;
121: PetscMalloc2(red->psubcomm->n*mlocal,PetscInt,&idx1,red->psubcomm->n*mlocal,PetscInt,&idx2);
122: j = 0;
123: for (k=0; k<red->psubcomm->n; k++) {
124: for (i=mstart; i<mend; i++) {
125: idx1[j] = i;
126: idx2[j++] = i + m*k;
127: }
128: }
129: ISCreateGeneral(comm,red->psubcomm->n*mlocal,idx1,PETSC_COPY_VALUES,&is1);
130: ISCreateGeneral(comm,red->psubcomm->n*mlocal,idx2,PETSC_COPY_VALUES,&is2);
131: VecScatterCreate(vec,is1,red->xdup,is2,&red->scatterin);
132: ISDestroy(&is1);
133: ISDestroy(&is2);
135: ISCreateStride(comm,mlocal,mstart+ red->psubcomm->color*m,1,&is1);
136: ISCreateStride(comm,mlocal,mstart,1,&is2);
137: VecScatterCreate(red->xdup,is1,vec,is2,&red->scatterout);
138: ISDestroy(&is1);
139: ISDestroy(&is2);
140: PetscFree2(idx1,idx2);
141: }
142: }
143: VecDestroy(&vec);
145: /* if pmatrix set by user is sequential then we do not need to gather the parallel matrix */
146: MPI_Comm_size(comm,&size);
147: if (size == 1) red->useparallelmat = PETSC_FALSE;
149: if (red->useparallelmat) {
150: if (pc->setupcalled == 1 && pc->flag == DIFFERENT_NONZERO_PATTERN) {
151: /* destroy old matrices */
152: MatDestroy(&red->pmats);
153: } else if (pc->setupcalled == 1) {
154: reuse = MAT_REUSE_MATRIX;
155: str = SAME_NONZERO_PATTERN;
156: }
158: /* grab the parallel matrix and put it into processors of a subcomminicator */
159: /*--------------------------------------------------------------------------*/
160: VecGetLocalSize(red->ysub,&mlocal_sub);
161: MatGetRedundantMatrix(pc->pmat,red->psubcomm->n,red->psubcomm->comm,mlocal_sub,reuse,&red->pmats);
162: /* tell PC of the subcommunicator its operators */
163: KSPSetOperators(red->ksp,red->pmats,red->pmats,str);
164: } else {
165: KSPSetOperators(red->ksp,pc->mat,pc->pmat,pc->flag);
166: }
167: if (pc->setfromoptionscalled) {
168: KSPSetFromOptions(red->ksp);
169: }
170: KSPSetUp(red->ksp);
171: return(0);
172: }
176: static PetscErrorCode PCApply_Redundant(PC pc,Vec x,Vec y)
177: {
178: PC_Redundant *red = (PC_Redundant*)pc->data;
180: PetscScalar *array;
183: /* scatter x to xdup */
184: VecScatterBegin(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD);
185: VecScatterEnd(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD);
187: /* place xdup's local array into xsub */
188: VecGetArray(red->xdup,&array);
189: VecPlaceArray(red->xsub,(const PetscScalar*)array);
191: /* apply preconditioner on each processor */
192: KSPSolve(red->ksp,red->xsub,red->ysub);
193: VecResetArray(red->xsub);
194: VecRestoreArray(red->xdup,&array);
196: /* place ysub's local array into ydup */
197: VecGetArray(red->ysub,&array);
198: VecPlaceArray(red->ydup,(const PetscScalar*)array);
200: /* scatter ydup to y */
201: VecScatterBegin(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD);
202: VecScatterEnd(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD);
203: VecResetArray(red->ydup);
204: VecRestoreArray(red->ysub,&array);
205: return(0);
206: }
210: static PetscErrorCode PCReset_Redundant(PC pc)
211: {
212: PC_Redundant *red = (PC_Redundant*)pc->data;
216: VecScatterDestroy(&red->scatterin);
217: VecScatterDestroy(&red->scatterout);
218: VecDestroy(&red->ysub);
219: VecDestroy(&red->xsub);
220: VecDestroy(&red->xdup);
221: VecDestroy(&red->ydup);
222: MatDestroy(&red->pmats);
223: if (red->ksp) {KSPReset(red->ksp);}
224: return(0);
225: }
229: static PetscErrorCode PCDestroy_Redundant(PC pc)
230: {
231: PC_Redundant *red = (PC_Redundant*)pc->data;
235: PCReset_Redundant(pc);
236: KSPDestroy(&red->ksp);
237: PetscSubcommDestroy(&red->psubcomm);
238: PetscFree(pc->data);
239: return(0);
240: }
244: static PetscErrorCode PCSetFromOptions_Redundant(PC pc)
245: {
247: PC_Redundant *red = (PC_Redundant*)pc->data;
250: PetscOptionsHead("Redundant options");
251: PetscOptionsInt("-pc_redundant_number","Number of redundant pc","PCRedundantSetNumber",red->nsubcomm,&red->nsubcomm,0);
252: PetscOptionsTail();
253: return(0);
254: }
258: static PetscErrorCode PCRedundantSetNumber_Redundant(PC pc,PetscInt nreds)
259: {
260: PC_Redundant *red = (PC_Redundant*)pc->data;
263: red->nsubcomm = nreds;
264: return(0);
265: }
269: /*@
270: PCRedundantSetNumber - Sets the number of redundant preconditioner contexts.
272: Logically Collective on PC
274: Input Parameters:
275: + pc - the preconditioner context
276: - nredundant - number of redundant preconditioner contexts; for example if you are using 64 MPI processes and
277: use an nredundant of 4 there will be 4 parallel solves each on 16 = 64/4 processes.
279: Level: advanced
281: .keywords: PC, redundant solve
282: @*/
283: PetscErrorCode PCRedundantSetNumber(PC pc,PetscInt nredundant)
284: {
289: if (nredundant <= 0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "num of redundant pc %D must be positive",nredundant);
290: PetscTryMethod(pc,"PCRedundantSetNumber_C",(PC,PetscInt),(pc,nredundant));
291: return(0);
292: }
296: static PetscErrorCode PCRedundantSetScatter_Redundant(PC pc,VecScatter in,VecScatter out)
297: {
298: PC_Redundant *red = (PC_Redundant*)pc->data;
302: PetscObjectReference((PetscObject)in);
303: VecScatterDestroy(&red->scatterin);
305: red->scatterin = in;
307: PetscObjectReference((PetscObject)out);
308: VecScatterDestroy(&red->scatterout);
309: red->scatterout = out;
310: return(0);
311: }
315: /*@
316: PCRedundantSetScatter - Sets the scatter used to copy values into the
317: redundant local solve and the scatter to move them back into the global
318: vector.
320: Logically Collective on PC
322: Input Parameters:
323: + pc - the preconditioner context
324: . in - the scatter to move the values in
325: - out - the scatter to move them out
327: Level: advanced
329: .keywords: PC, redundant solve
330: @*/
331: PetscErrorCode PCRedundantSetScatter(PC pc,VecScatter in,VecScatter out)
332: {
339: PetscTryMethod(pc,"PCRedundantSetScatter_C",(PC,VecScatter,VecScatter),(pc,in,out));
340: return(0);
341: }
345: static PetscErrorCode PCRedundantGetKSP_Redundant(PC pc,KSP *innerksp)
346: {
348: PC_Redundant *red = (PC_Redundant*)pc->data;
349: MPI_Comm comm,subcomm;
350: const char *prefix;
353: if (!red->psubcomm) {
354: PetscObjectGetComm((PetscObject)pc,&comm);
355: PetscSubcommCreate(comm,&red->psubcomm);
356: PetscSubcommSetNumber(red->psubcomm,red->nsubcomm);
357: PetscSubcommSetType(red->psubcomm,PETSC_SUBCOMM_INTERLACED);
358: PetscLogObjectMemory(pc,sizeof(PetscSubcomm));
360: /* create a new PC that processors in each subcomm have copy of */
361: subcomm = red->psubcomm->comm;
363: KSPCreate(subcomm,&red->ksp);
364: PetscObjectIncrementTabLevel((PetscObject)red->ksp,(PetscObject)pc,1);
365: PetscLogObjectParent(pc,red->ksp);
366: KSPSetType(red->ksp,KSPPREONLY);
367: KSPGetPC(red->ksp,&red->pc);
368: PCSetType(red->pc,PCLU);
370: PCGetOptionsPrefix(pc,&prefix);
371: KSPSetOptionsPrefix(red->ksp,prefix);
372: KSPAppendOptionsPrefix(red->ksp,"redundant_");
373: }
374: *innerksp = red->ksp;
375: return(0);
376: }
380: /*@
381: PCRedundantGetKSP - Gets the less parallel KSP created by the redundant PC.
383: Not Collective
385: Input Parameter:
386: . pc - the preconditioner context
388: Output Parameter:
389: . innerksp - the KSP on the smaller set of processes
391: Level: advanced
393: .keywords: PC, redundant solve
394: @*/
395: PetscErrorCode PCRedundantGetKSP(PC pc,KSP *innerksp)
396: {
402: PetscTryMethod(pc,"PCRedundantGetKSP_C",(PC,KSP*),(pc,innerksp));
403: return(0);
404: }
408: static PetscErrorCode PCRedundantGetOperators_Redundant(PC pc,Mat *mat,Mat *pmat)
409: {
410: PC_Redundant *red = (PC_Redundant*)pc->data;
413: if (mat) *mat = red->pmats;
414: if (pmat) *pmat = red->pmats;
415: return(0);
416: }
420: /*@
421: PCRedundantGetOperators - gets the sequential matrix and preconditioner matrix
423: Not Collective
425: Input Parameter:
426: . pc - the preconditioner context
428: Output Parameters:
429: + mat - the matrix
430: - pmat - the (possibly different) preconditioner matrix
432: Level: advanced
434: .keywords: PC, redundant solve
435: @*/
436: PetscErrorCode PCRedundantGetOperators(PC pc,Mat *mat,Mat *pmat)
437: {
444: PetscTryMethod(pc,"PCRedundantGetOperators_C",(PC,Mat*,Mat*),(pc,mat,pmat));
445: return(0);
446: }
448: /* -------------------------------------------------------------------------------------*/
449: /*MC
450: PCREDUNDANT - Runs a KSP solver with preconditioner for the entire problem on subgroups of processors
452: Options for the redundant preconditioners can be set with -redundant_pc_xxx for the redundant KSP with -redundant_ksp_xxx
454: Options Database:
455: . -pc_redundant_number <n> - number of redundant solves, for example if you are using 64 MPI processes and
456: use an n of 4 there will be 4 parallel solves each on 16 = 64/4 processes.
458: Level: intermediate
460: Notes: The default KSP is preonly and the default PC is LU.
462: Developer Notes: Note that PCSetInitialGuessNonzero() is not used by this class but likely should be.
464: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCRedundantSetScatter(),
465: PCRedundantGetKSP(), PCRedundantGetOperators(), PCRedundantSetNumber()
466: M*/
470: PETSC_EXTERN PetscErrorCode PCCreate_Redundant(PC pc)
471: {
473: PC_Redundant *red;
474: PetscMPIInt size;
477: PetscNewLog(pc,PC_Redundant,&red);
478: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
480: red->nsubcomm = size;
481: red->useparallelmat = PETSC_TRUE;
482: pc->data = (void*)red;
484: pc->ops->apply = PCApply_Redundant;
485: pc->ops->applytranspose = 0;
486: pc->ops->setup = PCSetUp_Redundant;
487: pc->ops->destroy = PCDestroy_Redundant;
488: pc->ops->reset = PCReset_Redundant;
489: pc->ops->setfromoptions = PCSetFromOptions_Redundant;
490: pc->ops->view = PCView_Redundant;
492: PetscObjectComposeFunction((PetscObject)pc,"PCRedundantSetScatter_C",PCRedundantSetScatter_Redundant);
493: PetscObjectComposeFunction((PetscObject)pc,"PCRedundantSetNumber_C",PCRedundantSetNumber_Redundant);
494: PetscObjectComposeFunction((PetscObject)pc,"PCRedundantGetKSP_C",PCRedundantGetKSP_Redundant);
495: PetscObjectComposeFunction((PetscObject)pc,"PCRedundantGetOperators_C",PCRedundantGetOperators_Redundant);
496: return(0);
497: }