Actual source code: redistribute.c

  1: /*
  2:   This file defines a "solve the problem redistributely on each subgroup of processor" preconditioner.
  3: */
  4: #include <petsc/private/pcimpl.h>
  5: #include <petscksp.h>

  7: typedef struct _PC_FieldSplitLink *PC_FieldSplitLink;
  8: struct _PC_FieldSplitLink {
  9:   char             *splitname;
 10:   IS                is;
 11:   PC_FieldSplitLink next, previous;
 12: };

 14: typedef struct {
 15:   KSP          ksp;
 16:   Vec          x, b;
 17:   VecScatter   scatter;
 18:   IS           is;
 19:   PetscInt     dcnt, *drows; /* these are the local rows that have only diagonal entry */
 20:   PetscScalar *diag;
 21:   Vec          work;
 22:   PetscBool    zerodiag;

 24:   PetscInt          nsplits;
 25:   PC_FieldSplitLink splitlinks;
 26: } PC_Redistribute;

 28: static PetscErrorCode PCFieldSplitSetIS_Redistribute(PC pc, const char splitname[], IS is)
 29: {
 30:   PC_Redistribute   *red  = (PC_Redistribute *)pc->data;
 31:   PC_FieldSplitLink *next = &red->splitlinks;

 33:   PetscFunctionBegin;
 34:   while (*next) next = &(*next)->next;
 35:   PetscCall(PetscNew(next));
 36:   if (splitname) {
 37:     PetscCall(PetscStrallocpy(splitname, &(*next)->splitname));
 38:   } else {
 39:     PetscCall(PetscMalloc1(8, &(*next)->splitname));
 40:     PetscCall(PetscSNPrintf((*next)->splitname, 7, "%" PetscInt_FMT, red->nsplits++));
 41:   }
 42:   PetscCall(PetscObjectReference((PetscObject)is));
 43:   PetscCall(ISDestroy(&(*next)->is));
 44:   (*next)->is = is;
 45:   PetscFunctionReturn(PETSC_SUCCESS);
 46: }

 48: static PetscErrorCode PCView_Redistribute(PC pc, PetscViewer viewer)
 49: {
 50:   PC_Redistribute *red = (PC_Redistribute *)pc->data;
 51:   PetscBool        iascii, isstring;
 52:   PetscInt         ncnt, N;

 54:   PetscFunctionBegin;
 55:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
 56:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
 57:   if (iascii) {
 58:     PetscCall(MPIU_Allreduce(&red->dcnt, &ncnt, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)pc)));
 59:     PetscCall(MatGetSize(pc->pmat, &N, NULL));
 60:     PetscCall(PetscViewerASCIIPrintf(viewer, "    Number rows eliminated %" PetscInt_FMT " Percentage rows eliminated %g\n", ncnt, (double)(100.0 * ((PetscReal)ncnt) / ((PetscReal)N))));
 61:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Redistribute preconditioner: \n"));
 62:     PetscCall(KSPView(red->ksp, viewer));
 63:   } else if (isstring) {
 64:     PetscCall(PetscViewerStringSPrintf(viewer, " Redistribute preconditioner"));
 65:     PetscCall(KSPView(red->ksp, viewer));
 66:   }
 67:   PetscFunctionReturn(PETSC_SUCCESS);
 68: }

 70: static PetscErrorCode PCSetUp_Redistribute(PC pc)
 71: {
 72:   PC_Redistribute         *red = (PC_Redistribute *)pc->data;
 73:   MPI_Comm                 comm;
 74:   PetscInt                 rstart, rend, nrstart, nrend, i, nz, cnt, *rows, ncnt, dcnt, *drows;
 75:   PetscLayout              map, nmap;
 76:   PetscMPIInt              size, tag, n;
 77:   PETSC_UNUSED PetscMPIInt imdex;
 78:   PetscInt                *source = NULL;
 79:   PetscMPIInt             *sizes  = NULL, nrecvs;
 80:   PetscInt                 j, nsends;
 81:   PetscInt                *owner = NULL, *starts = NULL, count, slen;
 82:   PetscInt                *rvalues, *svalues, recvtotal;
 83:   PetscMPIInt             *onodes1, *olengths1;
 84:   MPI_Request             *send_waits = NULL, *recv_waits = NULL;
 85:   MPI_Status               recv_status, *send_status;
 86:   Vec                      tvec, diag;
 87:   Mat                      tmat;
 88:   const PetscScalar       *d, *values;
 89:   const PetscInt          *cols;
 90:   PC_FieldSplitLink       *next = &red->splitlinks;

 92:   PetscFunctionBegin;
 93:   if (pc->setupcalled) {
 94:     PetscCheck(pc->flag == SAME_NONZERO_PATTERN, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "PC is not supported for a change in the nonzero structure of the matrix");
 95:     PetscCall(KSPGetOperators(red->ksp, NULL, &tmat));
 96:     PetscCall(MatCreateSubMatrix(pc->pmat, red->is, red->is, MAT_REUSE_MATRIX, &tmat));
 97:     PetscCall(KSPSetOperators(red->ksp, tmat, tmat));
 98:   } else {
 99:     PetscInt          NN;
100:     PC                ipc;
101:     PetscVoidFunction fptr;

103:     PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
104:     PetscCallMPI(MPI_Comm_size(comm, &size));
105:     PetscCall(PetscObjectGetNewTag((PetscObject)pc, &tag));

107:     /* count non-diagonal rows on process */
108:     PetscCall(MatGetOwnershipRange(pc->mat, &rstart, &rend));
109:     cnt = 0;
110:     for (i = rstart; i < rend; i++) {
111:       PetscCall(MatGetRow(pc->mat, i, &nz, &cols, &values));
112:       for (PetscInt j = 0; j < nz; j++) {
113:         if (values[j] != 0 && cols[j] != i) {
114:           cnt++;
115:           break;
116:         }
117:       }
118:       PetscCall(MatRestoreRow(pc->mat, i, &nz, &cols, &values));
119:     }
120:     PetscCall(PetscMalloc1(cnt, &rows));
121:     PetscCall(PetscMalloc1(rend - rstart - cnt, &drows));

123:     /* list non-diagonal rows on process */
124:     cnt  = 0;
125:     dcnt = 0;
126:     for (i = rstart; i < rend; i++) {
127:       PetscBool diagonly = PETSC_TRUE;
128:       PetscCall(MatGetRow(pc->mat, i, &nz, &cols, &values));
129:       for (PetscInt j = 0; j < nz; j++) {
130:         if (values[j] != 0 && cols[j] != i) {
131:           diagonly = PETSC_FALSE;
132:           break;
133:         }
134:       }
135:       if (!diagonly) rows[cnt++] = i;
136:       else drows[dcnt++] = i - rstart;
137:       PetscCall(MatRestoreRow(pc->mat, i, &nz, &cols, &values));
138:     }

140:     /* create PetscLayout for non-diagonal rows on each process */
141:     PetscCall(PetscLayoutCreate(comm, &map));
142:     PetscCall(PetscLayoutSetLocalSize(map, cnt));
143:     PetscCall(PetscLayoutSetBlockSize(map, 1));
144:     PetscCall(PetscLayoutSetUp(map));
145:     nrstart = map->rstart;
146:     nrend   = map->rend;

148:     /* create PetscLayout for load-balanced non-diagonal rows on each process */
149:     PetscCall(PetscLayoutCreate(comm, &nmap));
150:     PetscCall(MPIU_Allreduce(&cnt, &ncnt, 1, MPIU_INT, MPI_SUM, comm));
151:     PetscCall(PetscLayoutSetSize(nmap, ncnt));
152:     PetscCall(PetscLayoutSetBlockSize(nmap, 1));
153:     PetscCall(PetscLayoutSetUp(nmap));

155:     PetscCall(MatGetSize(pc->pmat, &NN, NULL));
156:     PetscCall(PetscInfo(pc, "Number of diagonal rows eliminated %" PetscInt_FMT ", percentage eliminated %g\n", NN - ncnt, (double)(((PetscReal)(NN - ncnt)) / ((PetscReal)(NN)))));

158:     if (size > 1) {
159:       /*
160:         the following block of code assumes MPI can send messages to self, which is not supported for MPI-uni hence we need to handle
161:         the size 1 case as a special case

163:        this code is taken from VecScatterCreate_PtoS()
164:        Determines what rows need to be moved where to
165:        load balance the non-diagonal rows
166:        */
167:       /*  count number of contributors to each processor */
168:       PetscCall(PetscMalloc2(size, &sizes, cnt, &owner));
169:       PetscCall(PetscArrayzero(sizes, size));
170:       j      = 0;
171:       nsends = 0;
172:       for (i = nrstart; i < nrend; i++) {
173:         if (i < nmap->range[j]) j = 0;
174:         for (; j < size; j++) {
175:           if (i < nmap->range[j + 1]) {
176:             if (!sizes[j]++) nsends++;
177:             owner[i - nrstart] = j;
178:             break;
179:           }
180:         }
181:       }
182:       /* inform other processors of number of messages and max length*/
183:       PetscCall(PetscGatherNumberOfMessages(comm, NULL, sizes, &nrecvs));
184:       PetscCall(PetscGatherMessageLengths(comm, nsends, nrecvs, sizes, &onodes1, &olengths1));
185:       PetscCall(PetscSortMPIIntWithArray(nrecvs, onodes1, olengths1));
186:       recvtotal = 0;
187:       for (i = 0; i < nrecvs; i++) recvtotal += olengths1[i];

189:       /* post receives:  rvalues - rows I will own; count - nu */
190:       PetscCall(PetscMalloc3(recvtotal, &rvalues, nrecvs, &source, nrecvs, &recv_waits));
191:       count = 0;
192:       for (i = 0; i < nrecvs; i++) {
193:         PetscCallMPI(MPI_Irecv((rvalues + count), olengths1[i], MPIU_INT, onodes1[i], tag, comm, recv_waits + i));
194:         count += olengths1[i];
195:       }

197:       /* do sends:
198:        1) starts[i] gives the starting index in svalues for stuff going to
199:        the ith processor
200:        */
201:       PetscCall(PetscMalloc3(cnt, &svalues, nsends, &send_waits, size, &starts));
202:       starts[0] = 0;
203:       for (i = 1; i < size; i++) starts[i] = starts[i - 1] + sizes[i - 1];
204:       for (i = 0; i < cnt; i++) svalues[starts[owner[i]]++] = rows[i];
205:       for (i = 0; i < cnt; i++) rows[i] = rows[i] - nrstart;
206:       red->drows = drows;
207:       red->dcnt  = dcnt;
208:       PetscCall(PetscFree(rows));

210:       starts[0] = 0;
211:       for (i = 1; i < size; i++) starts[i] = starts[i - 1] + sizes[i - 1];
212:       count = 0;
213:       for (i = 0; i < size; i++) {
214:         if (sizes[i]) PetscCallMPI(MPI_Isend(svalues + starts[i], sizes[i], MPIU_INT, i, tag, comm, send_waits + count++));
215:       }

217:       /*  wait on receives */
218:       count = nrecvs;
219:       slen  = 0;
220:       while (count) {
221:         PetscCallMPI(MPI_Waitany(nrecvs, recv_waits, &imdex, &recv_status));
222:         /* unpack receives into our local space */
223:         PetscCallMPI(MPI_Get_count(&recv_status, MPIU_INT, &n));
224:         slen += n;
225:         count--;
226:       }
227:       PetscCheck(slen == recvtotal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Total message lengths %" PetscInt_FMT " not expected %" PetscInt_FMT, slen, recvtotal);
228:       PetscCall(ISCreateGeneral(comm, slen, rvalues, PETSC_COPY_VALUES, &red->is));

230:       /* free all work space */
231:       PetscCall(PetscFree(olengths1));
232:       PetscCall(PetscFree(onodes1));
233:       PetscCall(PetscFree3(rvalues, source, recv_waits));
234:       PetscCall(PetscFree2(sizes, owner));
235:       if (nsends) { /* wait on sends */
236:         PetscCall(PetscMalloc1(nsends, &send_status));
237:         PetscCallMPI(MPI_Waitall(nsends, send_waits, send_status));
238:         PetscCall(PetscFree(send_status));
239:       }
240:       PetscCall(PetscFree3(svalues, send_waits, starts));
241:     } else {
242:       PetscCall(ISCreateGeneral(comm, cnt, rows, PETSC_OWN_POINTER, &red->is));
243:       red->drows = drows;
244:       red->dcnt  = dcnt;
245:       slen       = cnt;
246:     }
247:     PetscCall(PetscLayoutDestroy(&map));

249:     PetscCall(VecCreateMPI(comm, slen, PETSC_DETERMINE, &red->b));
250:     PetscCall(VecDuplicate(red->b, &red->x));
251:     PetscCall(MatCreateVecs(pc->pmat, &tvec, NULL));
252:     PetscCall(VecScatterCreate(tvec, red->is, red->b, NULL, &red->scatter));

254:     /* Map the PCFIELDSPLIT fields to redistributed KSP */
255:     PetscCall(KSPGetPC(red->ksp, &ipc));
256:     PetscCall(PetscObjectQueryFunction((PetscObject)ipc, "PCFieldSplitSetIS_C", &fptr));
257:     if (fptr && *next) {
258:       PetscScalar       *atvec;
259:       const PetscScalar *ab;
260:       PetscInt           primes[] = {2, 3, 5, 7, 11, 13, 17, 19};
261:       PetscInt           cnt      = 0;

263:       PetscCheck(red->nsplits <= (PetscInt)PETSC_STATIC_ARRAY_LENGTH(primes), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No support for this many fields");
264:       PetscCall(VecSet(tvec, 1.0));
265:       PetscCall(VecGetArray(tvec, &atvec));

267:       while (*next) {
268:         const PetscInt *indices;
269:         PetscInt        n;

271:         PetscCall(ISGetIndices((*next)->is, &indices));
272:         PetscCall(ISGetLocalSize((*next)->is, &n));
273:         for (PetscInt i = 0; i < n; i++) atvec[indices[i] - rstart] *= primes[cnt];
274:         PetscCall(ISRestoreIndices((*next)->is, &indices));
275:         cnt++;
276:         next = &(*next)->next;
277:       }
278:       PetscCall(VecRestoreArray(tvec, &atvec));
279:       PetscCall(VecScatterBegin(red->scatter, tvec, red->b, INSERT_VALUES, SCATTER_FORWARD));
280:       PetscCall(VecScatterEnd(red->scatter, tvec, red->b, INSERT_VALUES, SCATTER_FORWARD));
281:       cnt = 0;
282:       PetscCall(VecGetArrayRead(red->b, &ab));
283:       next = &red->splitlinks;
284:       while (*next) {
285:         PetscInt  n = 0;
286:         PetscInt *indices;
287:         IS        ris;

289:         for (PetscInt i = 0; i < nmap->rend - nmap->rstart; i++) {
290:           if (!(((PetscInt)PetscRealPart(ab[i])) % primes[cnt])) n++;
291:         }
292:         PetscCall(PetscMalloc1(n, &indices));
293:         n = 0;
294:         for (PetscInt i = 0; i < nmap->rend - nmap->rstart; i++) {
295:           if (!(((PetscInt)PetscRealPart(ab[i])) % primes[cnt])) indices[n++] = i + nmap->rstart;
296:         }
297:         PetscCall(ISCreateGeneral(comm, n, indices, PETSC_OWN_POINTER, &ris));
298:         PetscCall(PCFieldSplitSetIS(ipc, (*next)->splitname, ris));

300:         PetscCall(ISDestroy(&ris));
301:         cnt++;
302:         next = &(*next)->next;
303:       }
304:       PetscCall(VecRestoreArrayRead(red->b, &ab));
305:     }
306:     PetscCall(VecDestroy(&tvec));
307:     PetscCall(MatCreateSubMatrix(pc->pmat, red->is, red->is, MAT_INITIAL_MATRIX, &tmat));
308:     PetscCall(KSPSetOperators(red->ksp, tmat, tmat));
309:     PetscCall(MatDestroy(&tmat));
310:     PetscCall(PetscLayoutDestroy(&nmap));
311:   }

313:   /* get diagonal portion of matrix */
314:   PetscCall(PetscFree(red->diag));
315:   PetscCall(PetscMalloc1(red->dcnt, &red->diag));
316:   PetscCall(MatCreateVecs(pc->pmat, &diag, NULL));
317:   PetscCall(MatGetDiagonal(pc->pmat, diag));
318:   PetscCall(VecGetArrayRead(diag, &d));
319:   for (i = 0; i < red->dcnt; i++) {
320:     if (d[red->drows[i]] != 0) red->diag[i] = 1.0 / d[red->drows[i]];
321:     else {
322:       red->zerodiag = PETSC_TRUE;
323:       red->diag[i]  = 0.0;
324:     }
325:   }
326:   PetscCall(VecRestoreArrayRead(diag, &d));
327:   PetscCall(VecDestroy(&diag));
328:   PetscCall(KSPSetUp(red->ksp));
329:   PetscFunctionReturn(PETSC_SUCCESS);
330: }

332: static PetscErrorCode PCApply_Redistribute(PC pc, Vec b, Vec x)
333: {
334:   PC_Redistribute   *red   = (PC_Redistribute *)pc->data;
335:   PetscInt           dcnt  = red->dcnt, i;
336:   const PetscInt    *drows = red->drows;
337:   PetscScalar       *xwork;
338:   const PetscScalar *bwork, *diag = red->diag;
339:   PetscBool          nonzero_guess;

341:   PetscFunctionBegin;
342:   if (!red->work) PetscCall(VecDuplicate(b, &red->work));
343:   PetscCall(KSPGetInitialGuessNonzero(red->ksp, &nonzero_guess));
344:   if (nonzero_guess) {
345:     PetscCall(VecScatterBegin(red->scatter, x, red->x, INSERT_VALUES, SCATTER_FORWARD));
346:     PetscCall(VecScatterEnd(red->scatter, x, red->x, INSERT_VALUES, SCATTER_FORWARD));
347:   }

349:   /* compute the rows of solution that have diagonal entries only */
350:   PetscCall(VecSet(x, 0.0)); /* x = diag(A)^{-1} b */
351:   PetscCall(VecGetArray(x, &xwork));
352:   PetscCall(VecGetArrayRead(b, &bwork));
353:   if (red->zerodiag) {
354:     for (i = 0; i < dcnt; i++) {
355:       if (diag[i] == 0.0 && bwork[drows[i]] != 0.0) {
356:         PetscCheck(!pc->erroriffailure, PETSC_COMM_SELF, PETSC_ERR_CONV_FAILED, "Linear system is inconsistent, zero matrix row but nonzero right hand side");
357:         PetscCall(PetscInfo(pc, "Linear system is inconsistent, zero matrix row but nonzero right hand side\n"));
358:         PetscCall(VecSetInf(x));
359:         pc->failedreasonrank = PC_INCONSISTENT_RHS;
360:       }
361:     }
362:   }
363:   for (i = 0; i < dcnt; i++) xwork[drows[i]] = diag[i] * bwork[drows[i]];
364:   PetscCall(PetscLogFlops(dcnt));
365:   PetscCall(VecRestoreArray(red->work, &xwork));
366:   PetscCall(VecRestoreArrayRead(b, &bwork));
367:   /* update the right hand side for the reduced system with diagonal rows (and corresponding columns) removed */
368:   PetscCall(MatMult(pc->pmat, x, red->work));
369:   PetscCall(VecAYPX(red->work, -1.0, b)); /* red->work = b - A x */

371:   PetscCall(VecScatterBegin(red->scatter, red->work, red->b, INSERT_VALUES, SCATTER_FORWARD));
372:   PetscCall(VecScatterEnd(red->scatter, red->work, red->b, INSERT_VALUES, SCATTER_FORWARD));
373:   PetscCall(KSPSolve(red->ksp, red->b, red->x));
374:   PetscCall(KSPCheckSolve(red->ksp, pc, red->x));
375:   PetscCall(VecScatterBegin(red->scatter, red->x, x, INSERT_VALUES, SCATTER_REVERSE));
376:   PetscCall(VecScatterEnd(red->scatter, red->x, x, INSERT_VALUES, SCATTER_REVERSE));
377:   PetscFunctionReturn(PETSC_SUCCESS);
378: }

380: static PetscErrorCode PCDestroy_Redistribute(PC pc)
381: {
382:   PC_Redistribute  *red  = (PC_Redistribute *)pc->data;
383:   PC_FieldSplitLink next = red->splitlinks;

385:   PetscFunctionBegin;
386:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetIS_C", NULL));

388:   while (next) {
389:     PC_FieldSplitLink ilink;
390:     PetscCall(PetscFree(next->splitname));
391:     PetscCall(ISDestroy(&next->is));
392:     ilink = next;
393:     next  = next->next;
394:     PetscCall(PetscFree(ilink));
395:   }
396:   PetscCall(VecScatterDestroy(&red->scatter));
397:   PetscCall(ISDestroy(&red->is));
398:   PetscCall(VecDestroy(&red->b));
399:   PetscCall(VecDestroy(&red->x));
400:   PetscCall(KSPDestroy(&red->ksp));
401:   PetscCall(VecDestroy(&red->work));
402:   PetscCall(PetscFree(red->drows));
403:   PetscCall(PetscFree(red->diag));
404:   PetscCall(PetscFree(pc->data));
405:   PetscFunctionReturn(PETSC_SUCCESS);
406: }

408: static PetscErrorCode PCSetFromOptions_Redistribute(PC pc, PetscOptionItems *PetscOptionsObject)
409: {
410:   PC_Redistribute *red = (PC_Redistribute *)pc->data;

412:   PetscFunctionBegin;
413:   PetscCall(KSPSetFromOptions(red->ksp));
414:   PetscFunctionReturn(PETSC_SUCCESS);
415: }

417: /*@
418:   PCRedistributeGetKSP - Gets the `KSP` created by the `PCREDISTRIBUTE`

420:   Not Collective

422:   Input Parameter:
423: . pc - the preconditioner context

425:   Output Parameter:
426: . innerksp - the inner `KSP`

428:   Level: advanced

430: .seealso: [](ch_ksp), `KSP`, `PCREDISTRIBUTE`
431: @*/
432: PetscErrorCode PCRedistributeGetKSP(PC pc, KSP *innerksp)
433: {
434:   PC_Redistribute *red = (PC_Redistribute *)pc->data;

436:   PetscFunctionBegin;
438:   PetscAssertPointer(innerksp, 2);
439:   *innerksp = red->ksp;
440:   PetscFunctionReturn(PETSC_SUCCESS);
441: }

443: /*MC
444:      PCREDISTRIBUTE - Redistributes a matrix for load balancing, removing the rows (and the corresponding columns) that only have a diagonal entry and then
445:      applies a `KSP` to that new smaller matrix

447:      Level: intermediate

449:      Notes:
450:      Options for the redistribute `KSP` and `PC` with the options database prefix `-redistribute_`

452:      Usually run this with `-ksp_type preonly`

454:      If you have used `MatZeroRows()` to eliminate (for example, Dirichlet) boundary conditions for a symmetric problem then you can use, for example, `-ksp_type preonly
455:      -pc_type redistribute -redistribute_ksp_type cg -redistribute_pc_type bjacobi -redistribute_sub_pc_type icc` to take advantage of the symmetry.

457:      Supports the function `PCFieldSplitSetIS()`; pass the appropriate reduced field indices to an inner `PCFIELDSPLIT`, set with, for example
458:      `-ksp_type preonly -pc_type redistribute -redistribute_pc_type fieldsplit`. Does not support the `PCFIELDSPLIT` options database keys.

460:      This does NOT call a partitioner to reorder rows to lower communication; the ordering of the rows in the original matrix and redistributed matrix is the same. Rows are moved
461:      between MPI processes inside the preconditioner to balance the number of rows on each process.

463:      The matrix block information is lost with the possible removal of individual rows and columns of the matrix, thus the behavior of the preconditioner on the reduced
464:      system may be very different (worse) than running that preconditioner on the full system. This is specifically true for elasticity problems.

466:      Developer Note:
467:      Should add an option to this preconditioner to use a partitioner to redistribute the rows to lower communication.

469: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PCRedistributeGetKSP()`, `MatZeroRows()`, `PCFieldSplitSetIS()`, `PCFIELDSPLIT`
470: M*/

472: PETSC_EXTERN PetscErrorCode PCCreate_Redistribute(PC pc)
473: {
474:   PC_Redistribute *red;
475:   const char      *prefix;

477:   PetscFunctionBegin;
478:   PetscCall(PetscNew(&red));
479:   pc->data = (void *)red;

481:   pc->ops->apply          = PCApply_Redistribute;
482:   pc->ops->applytranspose = NULL;
483:   pc->ops->setup          = PCSetUp_Redistribute;
484:   pc->ops->destroy        = PCDestroy_Redistribute;
485:   pc->ops->setfromoptions = PCSetFromOptions_Redistribute;
486:   pc->ops->view           = PCView_Redistribute;

488:   PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &red->ksp));
489:   PetscCall(KSPSetNestLevel(red->ksp, pc->kspnestlevel));
490:   PetscCall(KSPSetErrorIfNotConverged(red->ksp, pc->erroriffailure));
491:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)red->ksp, (PetscObject)pc, 1));
492:   PetscCall(PCGetOptionsPrefix(pc, &prefix));
493:   PetscCall(KSPSetOptionsPrefix(red->ksp, prefix));
494:   PetscCall(KSPAppendOptionsPrefix(red->ksp, "redistribute_"));
495:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetIS_C", PCFieldSplitSetIS_Redistribute));
496:   PetscFunctionReturn(PETSC_SUCCESS);
497: }