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: }