Actual source code: redundant.c
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>
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: PetscBool shifttypeset;
19: MatFactorShiftType shifttype;
20: } PC_Redundant;
22: PetscErrorCode PCFactorSetShiftType_Redundant(PC pc,MatFactorShiftType shifttype)
23: {
24: PC_Redundant *red = (PC_Redundant*)pc->data;
26: if (red->ksp) {
27: PC pc;
28: KSPGetPC(red->ksp,&pc);
29: PCFactorSetShiftType(pc,shifttype);
30: } else {
31: red->shifttypeset = PETSC_TRUE;
32: red->shifttype = shifttype;
33: }
34: return 0;
35: }
37: static PetscErrorCode PCView_Redundant(PC pc,PetscViewer viewer)
38: {
39: PC_Redundant *red = (PC_Redundant*)pc->data;
40: PetscBool iascii,isstring;
41: PetscViewer subviewer;
43: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
44: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
45: if (iascii) {
46: if (!red->psubcomm) {
47: PetscViewerASCIIPrintf(viewer," Not yet setup\n");
48: } else {
49: PetscViewerASCIIPrintf(viewer," First (color=0) of %D PCs follows\n",red->nsubcomm);
50: PetscViewerGetSubViewer(viewer,((PetscObject)red->pc)->comm,&subviewer);
51: if (!red->psubcomm->color) { /* only view first redundant pc */
52: PetscViewerASCIIPushTab(subviewer);
53: KSPView(red->ksp,subviewer);
54: PetscViewerASCIIPopTab(subviewer);
55: }
56: PetscViewerRestoreSubViewer(viewer,((PetscObject)red->pc)->comm,&subviewer);
57: }
58: } else if (isstring) {
59: PetscViewerStringSPrintf(viewer," Redundant solver preconditioner");
60: }
61: return 0;
62: }
64: #include <../src/mat/impls/aij/mpi/mpiaij.h>
65: static PetscErrorCode PCSetUp_Redundant(PC pc)
66: {
67: PC_Redundant *red = (PC_Redundant*)pc->data;
68: PetscInt mstart,mend,mlocal,M;
69: PetscMPIInt size;
70: MPI_Comm comm,subcomm;
71: Vec x;
73: PetscObjectGetComm((PetscObject)pc,&comm);
75: /* if pmatrix set by user is sequential then we do not need to gather the parallel matrix */
76: MPI_Comm_size(comm,&size);
77: if (size == 1) red->useparallelmat = PETSC_FALSE;
79: if (!pc->setupcalled) {
80: PetscInt mloc_sub;
81: if (!red->psubcomm) { /* create red->psubcomm, new ksp and pc over subcomm */
82: KSP ksp;
83: PCRedundantGetKSP(pc,&ksp);
84: }
85: subcomm = PetscSubcommChild(red->psubcomm);
87: if (red->useparallelmat) {
88: /* grab the parallel matrix and put it into processors of a subcomminicator */
89: MatCreateRedundantMatrix(pc->pmat,red->psubcomm->n,subcomm,MAT_INITIAL_MATRIX,&red->pmats);
91: MPI_Comm_size(subcomm,&size);
92: if (size > 1) {
93: PetscBool foundpack,issbaij;
94: PetscObjectTypeCompare((PetscObject)red->pmats,MATMPISBAIJ,&issbaij);
95: if (!issbaij) {
96: MatGetFactorAvailable(red->pmats,NULL,MAT_FACTOR_LU,&foundpack);
97: } else {
98: MatGetFactorAvailable(red->pmats,NULL,MAT_FACTOR_CHOLESKY,&foundpack);
99: }
100: if (!foundpack) { /* reset default ksp and pc */
101: KSPSetType(red->ksp,KSPGMRES);
102: PCSetType(red->pc,PCBJACOBI);
103: } else {
104: PCFactorSetMatSolverType(red->pc,NULL);
105: }
106: }
108: KSPSetOperators(red->ksp,red->pmats,red->pmats);
110: /* get working vectors xsub and ysub */
111: MatCreateVecs(red->pmats,&red->xsub,&red->ysub);
113: /* create working vectors xdup and ydup.
114: xdup concatenates all xsub's contigously to form a mpi vector over dupcomm (see PetscSubcommCreate_interlaced())
115: ydup concatenates all ysub and has empty local arrays because ysub's arrays will be place into it.
116: Note: we use communicator dupcomm, not PetscObjectComm((PetscObject)pc)! */
117: MatGetLocalSize(red->pmats,&mloc_sub,NULL);
118: VecCreateMPI(PetscSubcommContiguousParent(red->psubcomm),mloc_sub,PETSC_DECIDE,&red->xdup);
119: VecCreateMPIWithArray(PetscSubcommContiguousParent(red->psubcomm),1,mloc_sub,PETSC_DECIDE,NULL,&red->ydup);
121: /* create vecscatters */
122: if (!red->scatterin) { /* efficiency of scatterin is independent from psubcomm_type! */
123: IS is1,is2;
124: PetscInt *idx1,*idx2,i,j,k;
126: MatCreateVecs(pc->pmat,&x,NULL);
127: VecGetSize(x,&M);
128: VecGetOwnershipRange(x,&mstart,&mend);
129: mlocal = mend - mstart;
130: PetscMalloc2(red->psubcomm->n*mlocal,&idx1,red->psubcomm->n*mlocal,&idx2);
131: j = 0;
132: for (k=0; k<red->psubcomm->n; k++) {
133: for (i=mstart; i<mend; i++) {
134: idx1[j] = i;
135: idx2[j++] = i + M*k;
136: }
137: }
138: ISCreateGeneral(comm,red->psubcomm->n*mlocal,idx1,PETSC_COPY_VALUES,&is1);
139: ISCreateGeneral(comm,red->psubcomm->n*mlocal,idx2,PETSC_COPY_VALUES,&is2);
140: VecScatterCreate(x,is1,red->xdup,is2,&red->scatterin);
141: ISDestroy(&is1);
142: ISDestroy(&is2);
144: /* Impl below is good for PETSC_SUBCOMM_INTERLACED (no inter-process communication) and PETSC_SUBCOMM_CONTIGUOUS (communication within subcomm) */
145: ISCreateStride(comm,mlocal,mstart+ red->psubcomm->color*M,1,&is1);
146: ISCreateStride(comm,mlocal,mstart,1,&is2);
147: VecScatterCreate(red->xdup,is1,x,is2,&red->scatterout);
148: ISDestroy(&is1);
149: ISDestroy(&is2);
150: PetscFree2(idx1,idx2);
151: VecDestroy(&x);
152: }
153: } else { /* !red->useparallelmat */
154: KSPSetOperators(red->ksp,pc->mat,pc->pmat);
155: }
156: } else { /* pc->setupcalled */
157: if (red->useparallelmat) {
158: MatReuse reuse;
159: /* grab the parallel matrix and put it into processors of a subcomminicator */
160: /*--------------------------------------------------------------------------*/
161: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
162: /* destroy old matrices */
163: MatDestroy(&red->pmats);
164: reuse = MAT_INITIAL_MATRIX;
165: } else {
166: reuse = MAT_REUSE_MATRIX;
167: }
168: MatCreateRedundantMatrix(pc->pmat,red->psubcomm->n,PetscSubcommChild(red->psubcomm),reuse,&red->pmats);
169: KSPSetOperators(red->ksp,red->pmats,red->pmats);
170: } else { /* !red->useparallelmat */
171: KSPSetOperators(red->ksp,pc->mat,pc->pmat);
172: }
173: }
175: if (pc->setfromoptionscalled) {
176: KSPSetFromOptions(red->ksp);
177: }
178: KSPSetUp(red->ksp);
179: return 0;
180: }
182: static PetscErrorCode PCApply_Redundant(PC pc,Vec x,Vec y)
183: {
184: PC_Redundant *red = (PC_Redundant*)pc->data;
185: PetscScalar *array;
187: if (!red->useparallelmat) {
188: KSPSolve(red->ksp,x,y);
189: KSPCheckSolve(red->ksp,pc,y);
190: return 0;
191: }
193: /* scatter x to xdup */
194: VecScatterBegin(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD);
195: VecScatterEnd(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD);
197: /* place xdup's local array into xsub */
198: VecGetArray(red->xdup,&array);
199: VecPlaceArray(red->xsub,(const PetscScalar*)array);
201: /* apply preconditioner on each processor */
202: KSPSolve(red->ksp,red->xsub,red->ysub);
203: KSPCheckSolve(red->ksp,pc,red->ysub);
204: VecResetArray(red->xsub);
205: VecRestoreArray(red->xdup,&array);
207: /* place ysub's local array into ydup */
208: VecGetArray(red->ysub,&array);
209: VecPlaceArray(red->ydup,(const PetscScalar*)array);
211: /* scatter ydup to y */
212: VecScatterBegin(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD);
213: VecScatterEnd(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD);
214: VecResetArray(red->ydup);
215: VecRestoreArray(red->ysub,&array);
216: return 0;
217: }
219: static PetscErrorCode PCApplyTranspose_Redundant(PC pc,Vec x,Vec y)
220: {
221: PC_Redundant *red = (PC_Redundant*)pc->data;
222: PetscScalar *array;
224: if (!red->useparallelmat) {
225: KSPSolveTranspose(red->ksp,x,y);
226: KSPCheckSolve(red->ksp,pc,y);
227: return 0;
228: }
230: /* scatter x to xdup */
231: VecScatterBegin(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD);
232: VecScatterEnd(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD);
234: /* place xdup's local array into xsub */
235: VecGetArray(red->xdup,&array);
236: VecPlaceArray(red->xsub,(const PetscScalar*)array);
238: /* apply preconditioner on each processor */
239: KSPSolveTranspose(red->ksp,red->xsub,red->ysub);
240: KSPCheckSolve(red->ksp,pc,red->ysub);
241: VecResetArray(red->xsub);
242: VecRestoreArray(red->xdup,&array);
244: /* place ysub's local array into ydup */
245: VecGetArray(red->ysub,&array);
246: VecPlaceArray(red->ydup,(const PetscScalar*)array);
248: /* scatter ydup to y */
249: VecScatterBegin(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD);
250: VecScatterEnd(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD);
251: VecResetArray(red->ydup);
252: VecRestoreArray(red->ysub,&array);
253: return 0;
254: }
256: static PetscErrorCode PCReset_Redundant(PC pc)
257: {
258: PC_Redundant *red = (PC_Redundant*)pc->data;
260: if (red->useparallelmat) {
261: VecScatterDestroy(&red->scatterin);
262: VecScatterDestroy(&red->scatterout);
263: VecDestroy(&red->ysub);
264: VecDestroy(&red->xsub);
265: VecDestroy(&red->xdup);
266: VecDestroy(&red->ydup);
267: }
268: MatDestroy(&red->pmats);
269: KSPReset(red->ksp);
270: return 0;
271: }
273: static PetscErrorCode PCDestroy_Redundant(PC pc)
274: {
275: PC_Redundant *red = (PC_Redundant*)pc->data;
277: PCReset_Redundant(pc);
278: KSPDestroy(&red->ksp);
279: PetscSubcommDestroy(&red->psubcomm);
280: PetscFree(pc->data);
281: return 0;
282: }
284: static PetscErrorCode PCSetFromOptions_Redundant(PetscOptionItems *PetscOptionsObject,PC pc)
285: {
286: PC_Redundant *red = (PC_Redundant*)pc->data;
288: PetscOptionsHead(PetscOptionsObject,"Redundant options");
289: PetscOptionsInt("-pc_redundant_number","Number of redundant pc","PCRedundantSetNumber",red->nsubcomm,&red->nsubcomm,NULL);
290: PetscOptionsTail();
291: return 0;
292: }
294: static PetscErrorCode PCRedundantSetNumber_Redundant(PC pc,PetscInt nreds)
295: {
296: PC_Redundant *red = (PC_Redundant*)pc->data;
298: red->nsubcomm = nreds;
299: return 0;
300: }
302: /*@
303: PCRedundantSetNumber - Sets the number of redundant preconditioner contexts.
305: Logically Collective on PC
307: Input Parameters:
308: + pc - the preconditioner context
309: - nredundant - number of redundant preconditioner contexts; for example if you are using 64 MPI processes and
310: use an nredundant of 4 there will be 4 parallel solves each on 16 = 64/4 processes.
312: Level: advanced
314: @*/
315: PetscErrorCode PCRedundantSetNumber(PC pc,PetscInt nredundant)
316: {
319: PetscTryMethod(pc,"PCRedundantSetNumber_C",(PC,PetscInt),(pc,nredundant));
320: return 0;
321: }
323: static PetscErrorCode PCRedundantSetScatter_Redundant(PC pc,VecScatter in,VecScatter out)
324: {
325: PC_Redundant *red = (PC_Redundant*)pc->data;
327: PetscObjectReference((PetscObject)in);
328: VecScatterDestroy(&red->scatterin);
330: red->scatterin = in;
332: PetscObjectReference((PetscObject)out);
333: VecScatterDestroy(&red->scatterout);
334: red->scatterout = out;
335: return 0;
336: }
338: /*@
339: PCRedundantSetScatter - Sets the scatter used to copy values into the
340: redundant local solve and the scatter to move them back into the global
341: vector.
343: Logically Collective on PC
345: Input Parameters:
346: + pc - the preconditioner context
347: . in - the scatter to move the values in
348: - out - the scatter to move them out
350: Level: advanced
352: @*/
353: PetscErrorCode PCRedundantSetScatter(PC pc,VecScatter in,VecScatter out)
354: {
358: PetscTryMethod(pc,"PCRedundantSetScatter_C",(PC,VecScatter,VecScatter),(pc,in,out));
359: return 0;
360: }
362: static PetscErrorCode PCRedundantGetKSP_Redundant(PC pc,KSP *innerksp)
363: {
364: PC_Redundant *red = (PC_Redundant*)pc->data;
365: MPI_Comm comm,subcomm;
366: const char *prefix;
367: PetscBool issbaij;
369: if (!red->psubcomm) {
370: PCGetOptionsPrefix(pc,&prefix);
372: PetscObjectGetComm((PetscObject)pc,&comm);
373: PetscSubcommCreate(comm,&red->psubcomm);
374: PetscSubcommSetNumber(red->psubcomm,red->nsubcomm);
375: PetscSubcommSetType(red->psubcomm,PETSC_SUBCOMM_CONTIGUOUS);
377: PetscSubcommSetOptionsPrefix(red->psubcomm,prefix);
378: PetscSubcommSetFromOptions(red->psubcomm);
379: PetscLogObjectMemory((PetscObject)pc,sizeof(PetscSubcomm));
381: /* create a new PC that processors in each subcomm have copy of */
382: subcomm = PetscSubcommChild(red->psubcomm);
384: KSPCreate(subcomm,&red->ksp);
385: KSPSetErrorIfNotConverged(red->ksp,pc->erroriffailure);
386: PetscObjectIncrementTabLevel((PetscObject)red->ksp,(PetscObject)pc,1);
387: PetscLogObjectParent((PetscObject)pc,(PetscObject)red->ksp);
388: KSPSetType(red->ksp,KSPPREONLY);
389: KSPGetPC(red->ksp,&red->pc);
390: PetscObjectTypeCompare((PetscObject)pc->pmat,MATSEQSBAIJ,&issbaij);
391: if (!issbaij) {
392: PetscObjectTypeCompare((PetscObject)pc->pmat,MATMPISBAIJ,&issbaij);
393: }
394: if (!issbaij) {
395: PCSetType(red->pc,PCLU);
396: } else {
397: PCSetType(red->pc,PCCHOLESKY);
398: }
399: if (red->shifttypeset) {
400: PCFactorSetShiftType(red->pc,red->shifttype);
401: red->shifttypeset = PETSC_FALSE;
402: }
403: KSPSetOptionsPrefix(red->ksp,prefix);
404: KSPAppendOptionsPrefix(red->ksp,"redundant_");
405: }
406: *innerksp = red->ksp;
407: return 0;
408: }
410: /*@
411: PCRedundantGetKSP - Gets the less parallel KSP created by the redundant PC.
413: Not Collective
415: Input Parameter:
416: . pc - the preconditioner context
418: Output Parameter:
419: . innerksp - the KSP on the smaller set of processes
421: Level: advanced
423: @*/
424: PetscErrorCode PCRedundantGetKSP(PC pc,KSP *innerksp)
425: {
428: PetscUseMethod(pc,"PCRedundantGetKSP_C",(PC,KSP*),(pc,innerksp));
429: return 0;
430: }
432: static PetscErrorCode PCRedundantGetOperators_Redundant(PC pc,Mat *mat,Mat *pmat)
433: {
434: PC_Redundant *red = (PC_Redundant*)pc->data;
436: if (mat) *mat = red->pmats;
437: if (pmat) *pmat = red->pmats;
438: return 0;
439: }
441: /*@
442: PCRedundantGetOperators - gets the sequential matrix and preconditioner matrix
444: Not Collective
446: Input Parameter:
447: . pc - the preconditioner context
449: Output Parameters:
450: + mat - the matrix
451: - pmat - the (possibly different) preconditioner matrix
453: Level: advanced
455: @*/
456: PetscErrorCode PCRedundantGetOperators(PC pc,Mat *mat,Mat *pmat)
457: {
461: PetscUseMethod(pc,"PCRedundantGetOperators_C",(PC,Mat*,Mat*),(pc,mat,pmat));
462: return 0;
463: }
465: /* -------------------------------------------------------------------------------------*/
466: /*MC
467: PCREDUNDANT - Runs a KSP solver with preconditioner for the entire problem on subgroups of processors
469: Options for the redundant preconditioners can be set with -redundant_pc_xxx for the redundant KSP with -redundant_ksp_xxx
471: Options Database:
472: . -pc_redundant_number <n> - number of redundant solves, for example if you are using 64 MPI processes and
473: use an n of 4 there will be 4 parallel solves each on 16 = 64/4 processes.
475: Level: intermediate
477: Notes:
478: The default KSP is preonly and the default PC is LU or CHOLESKY if Pmat is of type MATSBAIJ.
480: PCFactorSetShiftType() applied to this PC will convey they shift type into the inner PC if it is factorization based.
482: Developer Notes:
483: Note that PCSetInitialGuessNonzero() is not used by this class but likely should be.
485: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCRedundantSetScatter(),
486: PCRedundantGetKSP(), PCRedundantGetOperators(), PCRedundantSetNumber()
487: M*/
489: PETSC_EXTERN PetscErrorCode PCCreate_Redundant(PC pc)
490: {
491: PC_Redundant *red;
492: PetscMPIInt size;
494: PetscNewLog(pc,&red);
495: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
497: red->nsubcomm = size;
498: red->useparallelmat = PETSC_TRUE;
499: pc->data = (void*)red;
501: pc->ops->apply = PCApply_Redundant;
502: pc->ops->applytranspose = PCApplyTranspose_Redundant;
503: pc->ops->setup = PCSetUp_Redundant;
504: pc->ops->destroy = PCDestroy_Redundant;
505: pc->ops->reset = PCReset_Redundant;
506: pc->ops->setfromoptions = PCSetFromOptions_Redundant;
507: pc->ops->view = PCView_Redundant;
509: PetscObjectComposeFunction((PetscObject)pc,"PCRedundantSetScatter_C",PCRedundantSetScatter_Redundant);
510: PetscObjectComposeFunction((PetscObject)pc,"PCRedundantSetNumber_C",PCRedundantSetNumber_Redundant);
511: PetscObjectComposeFunction((PetscObject)pc,"PCRedundantGetKSP_C",PCRedundantGetKSP_Redundant);
512: PetscObjectComposeFunction((PetscObject)pc,"PCRedundantGetOperators_C",PCRedundantGetOperators_Redundant);
513: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftType_C",PCFactorSetShiftType_Redundant);
514: return 0;
515: }