Actual source code: redundant.c
petsc-3.3-p7 2013-05-11
2: /*
3: This file defines a "solve the problem redundantly on each subgroup of processor" preconditioner.
4: */
5: #include <petsc-private/pcimpl.h> /*I "petscpc.h" I*/
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 ((PetscObject)pc)->comm */
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: } else {
48: SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for PC redundant",((PetscObject)viewer)->type_name);
49: }
50: return(0);
51: }
55: static PetscErrorCode PCSetUp_Redundant(PC pc)
56: {
57: PC_Redundant *red = (PC_Redundant*)pc->data;
59: PetscInt mstart,mend,mlocal,m,mlocal_sub,rstart_sub,rend_sub,mloc_sub;
60: PetscMPIInt size;
61: MatReuse reuse = MAT_INITIAL_MATRIX;
62: MatStructure str = DIFFERENT_NONZERO_PATTERN;
63: MPI_Comm comm = ((PetscObject)pc)->comm,subcomm;
64: Vec vec;
65: PetscMPIInt subsize,subrank;
66: const char *prefix;
67: const PetscInt *range;
70: MatGetVecs(pc->pmat,&vec,0);
71: VecGetSize(vec,&m);
73: if (!pc->setupcalled) {
74: if (!red->psubcomm) {
75: PetscSubcommCreate(comm,&red->psubcomm);
76: PetscSubcommSetNumber(red->psubcomm,red->nsubcomm);
77: PetscSubcommSetType(red->psubcomm,PETSC_SUBCOMM_INTERLACED);
78: PetscLogObjectMemory(pc,sizeof(PetscSubcomm));
80: /* create a new PC that processors in each subcomm have copy of */
81: 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 {
93: subcomm = red->psubcomm->comm;
94: }
96: /* create working vectors xsub/ysub and xdup/ydup */
97: VecGetLocalSize(vec,&mlocal);
98: VecGetOwnershipRange(vec,&mstart,&mend);
100: /* get local size of xsub/ysub */
101: MPI_Comm_size(subcomm,&subsize);
102: MPI_Comm_rank(subcomm,&subrank);
103: MatGetOwnershipRanges(pc->pmat,&range);
104: rstart_sub = range[red->psubcomm->n*subrank]; /* rstart in xsub/ysub */
105: if (subrank+1 < subsize){
106: rend_sub = range[red->psubcomm->n*(subrank+1)];
107: } else {
108: rend_sub = m;
109: }
110: mloc_sub = rend_sub - rstart_sub;
111: VecCreateMPI(subcomm,mloc_sub,PETSC_DECIDE,&red->ysub);
112: /* create xsub with empty local arrays, because xdup's arrays will be placed into it */
113: VecCreateMPIWithArray(subcomm,1,mloc_sub,PETSC_DECIDE,PETSC_NULL,&red->xsub);
115: /* create xdup and ydup. ydup has empty local arrays because ysub's arrays will be place into it.
116: Note: we use communicator dupcomm, not ((PetscObject)pc)->comm! */
117: VecCreateMPI(red->psubcomm->dupparent,mloc_sub,PETSC_DECIDE,&red->xdup);
118: VecCreateMPIWithArray(red->psubcomm->dupparent,1,mloc_sub,PETSC_DECIDE,PETSC_NULL,&red->ydup);
119:
120: /* create vec scatters */
121: if (!red->scatterin){
122: IS is1,is2;
123: PetscInt *idx1,*idx2,i,j,k;
125: PetscMalloc2(red->psubcomm->n*mlocal,PetscInt,&idx1,red->psubcomm->n*mlocal,PetscInt,&idx2);
126: j = 0;
127: for (k=0; k<red->psubcomm->n; k++){
128: for (i=mstart; i<mend; i++){
129: idx1[j] = i;
130: idx2[j++] = i + m*k;
131: }
132: }
133: ISCreateGeneral(comm,red->psubcomm->n*mlocal,idx1,PETSC_COPY_VALUES,&is1);
134: ISCreateGeneral(comm,red->psubcomm->n*mlocal,idx2,PETSC_COPY_VALUES,&is2);
135: VecScatterCreate(vec,is1,red->xdup,is2,&red->scatterin);
136: ISDestroy(&is1);
137: ISDestroy(&is2);
139: ISCreateStride(comm,mlocal,mstart+ red->psubcomm->color*m,1,&is1);
140: ISCreateStride(comm,mlocal,mstart,1,&is2);
141: VecScatterCreate(red->xdup,is1,vec,is2,&red->scatterout);
142: ISDestroy(&is1);
143: ISDestroy(&is2);
144: PetscFree2(idx1,idx2);
145: }
146: }
147: VecDestroy(&vec);
149: /* if pmatrix set by user is sequential then we do not need to gather the parallel matrix */
150: MPI_Comm_size(comm,&size);
151: if (size == 1) {
152: red->useparallelmat = PETSC_FALSE;
153: }
155: if (red->useparallelmat) {
156: if (pc->setupcalled == 1 && pc->flag == DIFFERENT_NONZERO_PATTERN) {
157: /* destroy old matrices */
158: MatDestroy(&red->pmats);
159: } else if (pc->setupcalled == 1) {
160: reuse = MAT_REUSE_MATRIX;
161: str = SAME_NONZERO_PATTERN;
162: }
163:
164: /* grab the parallel matrix and put it into processors of a subcomminicator */
165: /*--------------------------------------------------------------------------*/
166: VecGetLocalSize(red->ysub,&mlocal_sub);
167: MatGetRedundantMatrix(pc->pmat,red->psubcomm->n,red->psubcomm->comm,mlocal_sub,reuse,&red->pmats);
168: /* tell PC of the subcommunicator its operators */
169: KSPSetOperators(red->ksp,red->pmats,red->pmats,str);
170: } else {
171: KSPSetOperators(red->ksp,pc->mat,pc->pmat,pc->flag);
172: }
173: if (pc->setfromoptionscalled){
174: KSPSetFromOptions(red->ksp);
175: }
176: KSPSetUp(red->ksp);
177: return(0);
178: }
182: static PetscErrorCode PCApply_Redundant(PC pc,Vec x,Vec y)
183: {
184: PC_Redundant *red = (PC_Redundant*)pc->data;
186: PetscScalar *array;
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);
192:
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);
201:
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 PCReset_Redundant(PC pc)
217: {
218: PC_Redundant *red = (PC_Redundant*)pc->data;
222: VecScatterDestroy(&red->scatterin);
223: VecScatterDestroy(&red->scatterout);
224: VecDestroy(&red->ysub);
225: VecDestroy(&red->xsub);
226: VecDestroy(&red->xdup);
227: VecDestroy(&red->ydup);
228: MatDestroy(&red->pmats);
229: if (red->ksp) {KSPReset(red->ksp);}
230: return(0);
231: }
235: static PetscErrorCode PCDestroy_Redundant(PC pc)
236: {
237: PC_Redundant *red = (PC_Redundant*)pc->data;
241: PCReset_Redundant(pc);
242: KSPDestroy(&red->ksp);
243: PetscSubcommDestroy(&red->psubcomm);
244: PetscFree(pc->data);
245: return(0);
246: }
250: static PetscErrorCode PCSetFromOptions_Redundant(PC pc)
251: {
253: PC_Redundant *red = (PC_Redundant*)pc->data;
256: PetscOptionsHead("Redundant options");
257: PetscOptionsInt("-pc_redundant_number","Number of redundant pc","PCRedundantSetNumber",red->nsubcomm,&red->nsubcomm,0);
258: PetscOptionsTail();
259: return(0);
260: }
262: EXTERN_C_BEGIN
265: PetscErrorCode PCRedundantSetNumber_Redundant(PC pc,PetscInt nreds)
266: {
267: PC_Redundant *red = (PC_Redundant*)pc->data;
270: red->nsubcomm = nreds;
271: return(0);
272: }
273: EXTERN_C_END
277: /*@
278: PCRedundantSetNumber - Sets the number of redundant preconditioner contexts.
280: Logically Collective on PC
282: Input Parameters:
283: + pc - the preconditioner context
284: - nredundant - number of redundant preconditioner contexts; for example if you are using 64 MPI processes and
285: use an nredundant of 4 there will be 4 parallel solves each on 16 = 64/4 processes.
287: Level: advanced
289: .keywords: PC, redundant solve
290: @*/
291: PetscErrorCode PCRedundantSetNumber(PC pc,PetscInt nredundant)
292: {
297: if (nredundant <= 0) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONG, "num of redundant pc %D must be positive",nredundant);
298: PetscTryMethod(pc,"PCRedundantSetNumber_C",(PC,PetscInt),(pc,nredundant));
299: return(0);
300: }
302: EXTERN_C_BEGIN
305: PetscErrorCode PCRedundantSetScatter_Redundant(PC pc,VecScatter in,VecScatter out)
306: {
307: PC_Redundant *red = (PC_Redundant*)pc->data;
311: PetscObjectReference((PetscObject)in);
312: VecScatterDestroy(&red->scatterin);
313: red->scatterin = in;
314: PetscObjectReference((PetscObject)out);
315: VecScatterDestroy(&red->scatterout);
316: red->scatterout = out;
317: return(0);
318: }
319: EXTERN_C_END
323: /*@
324: PCRedundantSetScatter - Sets the scatter used to copy values into the
325: redundant local solve and the scatter to move them back into the global
326: vector.
328: Logically Collective on PC
330: Input Parameters:
331: + pc - the preconditioner context
332: . in - the scatter to move the values in
333: - out - the scatter to move them out
335: Level: advanced
337: .keywords: PC, redundant solve
338: @*/
339: PetscErrorCode PCRedundantSetScatter(PC pc,VecScatter in,VecScatter out)
340: {
347: PetscTryMethod(pc,"PCRedundantSetScatter_C",(PC,VecScatter,VecScatter),(pc,in,out));
348: return(0);
349: }
351: EXTERN_C_BEGIN
354: PetscErrorCode PCRedundantGetKSP_Redundant(PC pc,KSP *innerksp)
355: {
357: PC_Redundant *red = (PC_Redundant*)pc->data;
358: MPI_Comm comm,subcomm;
359: const char *prefix;
362: if (!red->psubcomm) {
363: PetscObjectGetComm((PetscObject)pc,&comm);
364: PetscSubcommCreate(comm,&red->psubcomm);
365: PetscSubcommSetNumber(red->psubcomm,red->nsubcomm);
366: PetscSubcommSetType(red->psubcomm,PETSC_SUBCOMM_INTERLACED);
367: PetscLogObjectMemory(pc,sizeof(PetscSubcomm));
369: /* create a new PC that processors in each subcomm have copy of */
370: subcomm = red->psubcomm->comm;
371: KSPCreate(subcomm,&red->ksp);
372: PetscObjectIncrementTabLevel((PetscObject)red->ksp,(PetscObject)pc,1);
373: PetscLogObjectParent(pc,red->ksp);
374: KSPSetType(red->ksp,KSPPREONLY);
375: KSPGetPC(red->ksp,&red->pc);
376: PCSetType(red->pc,PCLU);
378: PCGetOptionsPrefix(pc,&prefix);
379: KSPSetOptionsPrefix(red->ksp,prefix);
380: KSPAppendOptionsPrefix(red->ksp,"redundant_");
381: }
382: *innerksp = red->ksp;
383: return(0);
384: }
385: EXTERN_C_END
389: /*@
390: PCRedundantGetKSP - Gets the less parallel KSP created by the redundant PC.
392: Not Collective
394: Input Parameter:
395: . pc - the preconditioner context
397: Output Parameter:
398: . innerksp - the KSP on the smaller set of processes
400: Level: advanced
402: .keywords: PC, redundant solve
403: @*/
404: PetscErrorCode PCRedundantGetKSP(PC pc,KSP *innerksp)
405: {
411: PetscTryMethod(pc,"PCRedundantGetKSP_C",(PC,KSP*),(pc,innerksp));
412: return(0);
413: }
415: EXTERN_C_BEGIN
418: PetscErrorCode PCRedundantGetOperators_Redundant(PC pc,Mat *mat,Mat *pmat)
419: {
420: PC_Redundant *red = (PC_Redundant*)pc->data;
423: if (mat) *mat = red->pmats;
424: if (pmat) *pmat = red->pmats;
425: return(0);
426: }
427: EXTERN_C_END
431: /*@
432: PCRedundantGetOperators - gets the sequential matrix and preconditioner matrix
434: Not Collective
436: Input Parameter:
437: . pc - the preconditioner context
439: Output Parameters:
440: + mat - the matrix
441: - pmat - the (possibly different) preconditioner matrix
443: Level: advanced
445: .keywords: PC, redundant solve
446: @*/
447: PetscErrorCode PCRedundantGetOperators(PC pc,Mat *mat,Mat *pmat)
448: {
455: PetscTryMethod(pc,"PCRedundantGetOperators_C",(PC,Mat*,Mat*),(pc,mat,pmat));
456: return(0);
457: }
459: /* -------------------------------------------------------------------------------------*/
460: /*MC
461: PCREDUNDANT - Runs a KSP solver with preconditioner for the entire problem on subgroups of processors
463: Options for the redundant preconditioners can be set with -redundant_pc_xxx for the redundant KSP with -redundant_ksp_xxx
465: Options Database:
466: . -pc_redundant_number <n> - number of redundant solves, for example if you are using 64 MPI processes and
467: use an n of 4 there will be 4 parallel solves each on 16 = 64/4 processes.
469: Level: intermediate
471: Notes: The default KSP is preonly and the default PC is LU.
473: Developer Notes: Note that PCSetInitialGuessNonzero() is not used by this class but likely should be.
475: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCRedundantSetScatter(),
476: PCRedundantGetKSP(), PCRedundantGetOperators(), PCRedundantSetNumber()
477: M*/
479: EXTERN_C_BEGIN
482: PetscErrorCode PCCreate_Redundant(PC pc)
483: {
485: PC_Redundant *red;
486: PetscMPIInt size;
487:
489: PetscNewLog(pc,PC_Redundant,&red);
490: MPI_Comm_size(((PetscObject)pc)->comm,&size);
491: red->nsubcomm = size;
492: red->useparallelmat = PETSC_TRUE;
493: pc->data = (void*)red;
495: pc->ops->apply = PCApply_Redundant;
496: pc->ops->applytranspose = 0;
497: pc->ops->setup = PCSetUp_Redundant;
498: pc->ops->destroy = PCDestroy_Redundant;
499: pc->ops->reset = PCReset_Redundant;
500: pc->ops->setfromoptions = PCSetFromOptions_Redundant;
501: pc->ops->view = PCView_Redundant;
502: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCRedundantSetScatter_C","PCRedundantSetScatter_Redundant",
503: PCRedundantSetScatter_Redundant);
504: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCRedundantSetNumber_C","PCRedundantSetNumber_Redundant",
505: PCRedundantSetNumber_Redundant);
506: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCRedundantGetKSP_C","PCRedundantGetKSP_Redundant",
507: PCRedundantGetKSP_Redundant);
508: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCRedundantGetOperators_C","PCRedundantGetOperators_Redundant",
509: PCRedundantGetOperators_Redundant);
510: return(0);
511: }
512: EXTERN_C_END