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 %" PetscInt_FMT " 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) KSPSetFromOptions(red->ksp);
176: KSPSetUp(red->ksp);
177: return 0;
178: }
180: static PetscErrorCode PCApply_Redundant(PC pc, Vec x, Vec y)
181: {
182: PC_Redundant *red = (PC_Redundant *)pc->data;
183: PetscScalar *array;
185: if (!red->useparallelmat) {
186: KSPSolve(red->ksp, x, y);
187: KSPCheckSolve(red->ksp, pc, y);
188: return 0;
189: }
191: /* scatter x to xdup */
192: VecScatterBegin(red->scatterin, x, red->xdup, INSERT_VALUES, SCATTER_FORWARD);
193: VecScatterEnd(red->scatterin, x, red->xdup, INSERT_VALUES, SCATTER_FORWARD);
195: /* place xdup's local array into xsub */
196: VecGetArray(red->xdup, &array);
197: VecPlaceArray(red->xsub, (const PetscScalar *)array);
199: /* apply preconditioner on each processor */
200: KSPSolve(red->ksp, red->xsub, red->ysub);
201: KSPCheckSolve(red->ksp, pc, red->ysub);
202: VecResetArray(red->xsub);
203: VecRestoreArray(red->xdup, &array);
205: /* place ysub's local array into ydup */
206: VecGetArray(red->ysub, &array);
207: VecPlaceArray(red->ydup, (const PetscScalar *)array);
209: /* scatter ydup to y */
210: VecScatterBegin(red->scatterout, red->ydup, y, INSERT_VALUES, SCATTER_FORWARD);
211: VecScatterEnd(red->scatterout, red->ydup, y, INSERT_VALUES, SCATTER_FORWARD);
212: VecResetArray(red->ydup);
213: VecRestoreArray(red->ysub, &array);
214: return 0;
215: }
217: static PetscErrorCode PCApplyTranspose_Redundant(PC pc, Vec x, Vec y)
218: {
219: PC_Redundant *red = (PC_Redundant *)pc->data;
220: PetscScalar *array;
222: if (!red->useparallelmat) {
223: KSPSolveTranspose(red->ksp, x, y);
224: KSPCheckSolve(red->ksp, pc, 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: KSPCheckSolve(red->ksp, pc, red->ysub);
239: VecResetArray(red->xsub);
240: VecRestoreArray(red->xdup, &array);
242: /* place ysub's local array into ydup */
243: VecGetArray(red->ysub, &array);
244: VecPlaceArray(red->ydup, (const PetscScalar *)array);
246: /* scatter ydup to y */
247: VecScatterBegin(red->scatterout, red->ydup, y, INSERT_VALUES, SCATTER_FORWARD);
248: VecScatterEnd(red->scatterout, red->ydup, y, INSERT_VALUES, SCATTER_FORWARD);
249: VecResetArray(red->ydup);
250: VecRestoreArray(red->ysub, &array);
251: return 0;
252: }
254: static PetscErrorCode PCReset_Redundant(PC pc)
255: {
256: PC_Redundant *red = (PC_Redundant *)pc->data;
258: if (red->useparallelmat) {
259: VecScatterDestroy(&red->scatterin);
260: VecScatterDestroy(&red->scatterout);
261: VecDestroy(&red->ysub);
262: VecDestroy(&red->xsub);
263: VecDestroy(&red->xdup);
264: VecDestroy(&red->ydup);
265: }
266: MatDestroy(&red->pmats);
267: KSPReset(red->ksp);
268: return 0;
269: }
271: static PetscErrorCode PCDestroy_Redundant(PC pc)
272: {
273: PC_Redundant *red = (PC_Redundant *)pc->data;
275: PCReset_Redundant(pc);
276: KSPDestroy(&red->ksp);
277: PetscSubcommDestroy(&red->psubcomm);
278: PetscObjectComposeFunction((PetscObject)pc, "PCRedundantSetScatter_C", NULL);
279: PetscObjectComposeFunction((PetscObject)pc, "PCRedundantSetNumber_C", NULL);
280: PetscObjectComposeFunction((PetscObject)pc, "PCRedundantGetKSP_C", NULL);
281: PetscObjectComposeFunction((PetscObject)pc, "PCRedundantGetOperators_C", NULL);
282: PetscObjectComposeFunction((PetscObject)pc, "PCFactorSetShiftType_C", NULL);
283: PetscFree(pc->data);
284: return 0;
285: }
287: static PetscErrorCode PCSetFromOptions_Redundant(PC pc, PetscOptionItems *PetscOptionsObject)
288: {
289: PC_Redundant *red = (PC_Redundant *)pc->data;
291: PetscOptionsHeadBegin(PetscOptionsObject, "Redundant options");
292: PetscOptionsInt("-pc_redundant_number", "Number of redundant pc", "PCRedundantSetNumber", red->nsubcomm, &red->nsubcomm, NULL);
293: PetscOptionsHeadEnd();
294: return 0;
295: }
297: static PetscErrorCode PCRedundantSetNumber_Redundant(PC pc, PetscInt nreds)
298: {
299: PC_Redundant *red = (PC_Redundant *)pc->data;
301: red->nsubcomm = nreds;
302: return 0;
303: }
305: /*@
306: PCRedundantSetNumber - Sets the number of redundant preconditioner contexts.
308: Logically Collective
310: Input Parameters:
311: + pc - the preconditioner context
312: - nredundant - number of redundant preconditioner contexts; for example if you are using 64 MPI processes and
313: use an nredundant of 4 there will be 4 parallel solves each on 16 = 64/4 processes.
315: Level: advanced
317: .seealso: `PCREDUNDANT`
318: @*/
319: PetscErrorCode PCRedundantSetNumber(PC pc, PetscInt nredundant)
320: {
323: PetscTryMethod(pc, "PCRedundantSetNumber_C", (PC, PetscInt), (pc, nredundant));
324: return 0;
325: }
327: static PetscErrorCode PCRedundantSetScatter_Redundant(PC pc, VecScatter in, VecScatter out)
328: {
329: PC_Redundant *red = (PC_Redundant *)pc->data;
331: PetscObjectReference((PetscObject)in);
332: VecScatterDestroy(&red->scatterin);
334: red->scatterin = in;
336: PetscObjectReference((PetscObject)out);
337: VecScatterDestroy(&red->scatterout);
338: red->scatterout = out;
339: return 0;
340: }
342: /*@
343: PCRedundantSetScatter - Sets the scatter used to copy values into the
344: redundant local solve and the scatter to move them back into the global
345: vector.
347: Logically Collective
349: Input Parameters:
350: + pc - the preconditioner context
351: . in - the scatter to move the values in
352: - out - the scatter to move them out
354: Level: advanced
356: .seealso: `PCREDUNDANT`
357: @*/
358: PetscErrorCode PCRedundantSetScatter(PC pc, VecScatter in, VecScatter out)
359: {
363: PetscTryMethod(pc, "PCRedundantSetScatter_C", (PC, VecScatter, VecScatter), (pc, in, out));
364: return 0;
365: }
367: static PetscErrorCode PCRedundantGetKSP_Redundant(PC pc, KSP *innerksp)
368: {
369: PC_Redundant *red = (PC_Redundant *)pc->data;
370: MPI_Comm comm, subcomm;
371: const char *prefix;
372: PetscBool issbaij;
374: if (!red->psubcomm) {
375: PCGetOptionsPrefix(pc, &prefix);
377: PetscObjectGetComm((PetscObject)pc, &comm);
378: PetscSubcommCreate(comm, &red->psubcomm);
379: PetscSubcommSetNumber(red->psubcomm, red->nsubcomm);
380: PetscSubcommSetType(red->psubcomm, PETSC_SUBCOMM_CONTIGUOUS);
382: PetscSubcommSetOptionsPrefix(red->psubcomm, prefix);
383: PetscSubcommSetFromOptions(red->psubcomm);
385: /* create a new PC that processors in each subcomm have copy of */
386: subcomm = PetscSubcommChild(red->psubcomm);
388: KSPCreate(subcomm, &red->ksp);
389: KSPSetErrorIfNotConverged(red->ksp, pc->erroriffailure);
390: PetscObjectIncrementTabLevel((PetscObject)red->ksp, (PetscObject)pc, 1);
391: KSPSetType(red->ksp, KSPPREONLY);
392: KSPGetPC(red->ksp, &red->pc);
393: PetscObjectTypeCompare((PetscObject)pc->pmat, MATSEQSBAIJ, &issbaij);
394: if (!issbaij) PetscObjectTypeCompare((PetscObject)pc->pmat, MATMPISBAIJ, &issbaij);
395: if (!issbaij) {
396: PCSetType(red->pc, PCLU);
397: } else {
398: PCSetType(red->pc, PCCHOLESKY);
399: }
400: if (red->shifttypeset) {
401: PCFactorSetShiftType(red->pc, red->shifttype);
402: red->shifttypeset = PETSC_FALSE;
403: }
404: KSPSetOptionsPrefix(red->ksp, prefix);
405: KSPAppendOptionsPrefix(red->ksp, "redundant_");
406: }
407: *innerksp = red->ksp;
408: return 0;
409: }
411: /*@
412: PCRedundantGetKSP - Gets the less parallel `KSP` created by the redundant `PC`.
414: Not Collective
416: Input Parameter:
417: . pc - the preconditioner context
419: Output Parameter:
420: . innerksp - the `KSP` on the smaller set of processes
422: Level: advanced
424: .seealso: `PCREDUNDANT`
425: @*/
426: PetscErrorCode PCRedundantGetKSP(PC pc, KSP *innerksp)
427: {
430: PetscUseMethod(pc, "PCRedundantGetKSP_C", (PC, KSP *), (pc, innerksp));
431: return 0;
432: }
434: static PetscErrorCode PCRedundantGetOperators_Redundant(PC pc, Mat *mat, Mat *pmat)
435: {
436: PC_Redundant *red = (PC_Redundant *)pc->data;
438: if (mat) *mat = red->pmats;
439: if (pmat) *pmat = red->pmats;
440: return 0;
441: }
443: /*@
444: PCRedundantGetOperators - gets the sequential matrix and preconditioner matrix
446: Not Collective
448: Input Parameter:
449: . pc - the preconditioner context
451: Output Parameters:
452: + mat - the matrix
453: - pmat - the (possibly different) preconditioner matrix
455: Level: advanced
457: .seealso: `PCREDUNDANT`
458: @*/
459: PetscErrorCode PCRedundantGetOperators(PC pc, Mat *mat, Mat *pmat)
460: {
464: PetscUseMethod(pc, "PCRedundantGetOperators_C", (PC, Mat *, Mat *), (pc, mat, pmat));
465: return 0;
466: }
468: /*MC
469: PCREDUNDANT - Runs a `KSP` solver with preconditioner for the entire problem on subgroups of processors
471: Options Database Key:
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: Options for the redundant preconditioners can be set using the options database prefix -redundant_
480: The default `KSP` is preonly and the default `PC` is `PCLU` or `PCCHOLESKY` if Pmat is of type `MATSBAIJ`.
482: `PCFactorSetShiftType()` applied to this `PC` will convey they shift type into the inner `PC` if it is factorization based.
484: Developer Note:
485: `PCSetInitialGuessNonzero()` is not used by this class but likely should be.
487: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PCRedundantSetScatter()`,
488: `PCRedundantGetKSP()`, `PCRedundantGetOperators()`, `PCRedundantSetNumber()`, `PCREDISTRIBUTE`
489: M*/
491: PETSC_EXTERN PetscErrorCode PCCreate_Redundant(PC pc)
492: {
493: PC_Redundant *red;
494: PetscMPIInt size;
496: PetscNew(&red);
497: MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size);
499: red->nsubcomm = size;
500: red->useparallelmat = PETSC_TRUE;
501: pc->data = (void *)red;
503: pc->ops->apply = PCApply_Redundant;
504: pc->ops->applytranspose = PCApplyTranspose_Redundant;
505: pc->ops->setup = PCSetUp_Redundant;
506: pc->ops->destroy = PCDestroy_Redundant;
507: pc->ops->reset = PCReset_Redundant;
508: pc->ops->setfromoptions = PCSetFromOptions_Redundant;
509: pc->ops->view = PCView_Redundant;
511: PetscObjectComposeFunction((PetscObject)pc, "PCRedundantSetScatter_C", PCRedundantSetScatter_Redundant);
512: PetscObjectComposeFunction((PetscObject)pc, "PCRedundantSetNumber_C", PCRedundantSetNumber_Redundant);
513: PetscObjectComposeFunction((PetscObject)pc, "PCRedundantGetKSP_C", PCRedundantGetKSP_Redundant);
514: PetscObjectComposeFunction((PetscObject)pc, "PCRedundantGetOperators_C", PCRedundantGetOperators_Redundant);
515: PetscObjectComposeFunction((PetscObject)pc, "PCFactorSetShiftType_C", PCFactorSetShiftType_Redundant);
516: return 0;
517: }