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:     PetscCallMPI(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, 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, nsends;
     80:   PetscInt                 j;
     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:     PetscVoidFn *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 (PetscInt 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 (PetscInt 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:     PetscCallMPI(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 (PetscInt 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 (PetscMPIInt 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 (PetscMPIInt i = 0; i < nrecvs; i++) {
    193:         PetscCallMPI(MPIU_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 (PetscMPIInt i = 1; i < size; i++) starts[i] = starts[i - 1] + sizes[i - 1];
    204:       for (PetscInt i = 0; i < cnt; i++) svalues[starts[owner[i]]++] = rows[i];
    205:       for (PetscInt 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 (PetscMPIInt i = 1; i < size; i++) starts[i] = starts[i - 1] + sizes[i - 1];
    212:       count = 0;
    213:       for (PetscMPIInt i = 0; i < size; i++) {
    214:         if (sizes[i]) PetscCallMPI(MPIU_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 (PetscInt 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:         pc->failedreasonrank = PC_INCONSISTENT_RHS;
    359:       }
    360:     }
    361:     PetscCall(VecFlag(x, pc->failedreasonrank == PC_INCONSISTENT_RHS));
    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 PCApplyTranspose_Redistribute(PC pc, Vec b, Vec x)
    381: {
    382:   PC_Redistribute   *red   = (PC_Redistribute *)pc->data;
    383:   PetscInt           dcnt  = red->dcnt, i;
    384:   const PetscInt    *drows = red->drows;
    385:   PetscScalar       *xwork;
    386:   const PetscScalar *bwork, *diag = red->diag;
    387:   PetscBool          set, flg     = PETSC_FALSE, nonzero_guess;

    389:   PetscFunctionBegin;
    390:   PetscCall(MatIsStructurallySymmetricKnown(pc->pmat, &set, &flg));
    391:   PetscCheck(set || flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "PCApplyTranspose() not implemented for structurally unsymmetric Mat");
    392:   if (!red->work) PetscCall(VecDuplicate(b, &red->work));
    393:   PetscCall(KSPGetInitialGuessNonzero(red->ksp, &nonzero_guess));
    394:   if (nonzero_guess) {
    395:     PetscCall(VecScatterBegin(red->scatter, x, red->x, INSERT_VALUES, SCATTER_FORWARD));
    396:     PetscCall(VecScatterEnd(red->scatter, x, red->x, INSERT_VALUES, SCATTER_FORWARD));
    397:   }

    399:   /* compute the rows of solution that have diagonal entries only */
    400:   PetscCall(VecSet(x, 0.0)); /* x = diag(A)^{-1} b */
    401:   PetscCall(VecGetArray(x, &xwork));
    402:   PetscCall(VecGetArrayRead(b, &bwork));
    403:   if (red->zerodiag) {
    404:     for (i = 0; i < dcnt; i++) {
    405:       if (diag[i] == 0.0 && bwork[drows[i]] != 0.0) {
    406:         PetscCheck(!pc->erroriffailure, PETSC_COMM_SELF, PETSC_ERR_CONV_FAILED, "Linear system is inconsistent, zero matrix row but nonzero right-hand side");
    407:         PetscCall(PetscInfo(pc, "Linear system is inconsistent, zero matrix row but nonzero right-hand side\n"));
    408:         pc->failedreasonrank = PC_INCONSISTENT_RHS;
    409:       }
    410:     }
    411:     PetscCall(VecFlag(x, pc->failedreasonrank == PC_INCONSISTENT_RHS));
    412:   }
    413:   for (i = 0; i < dcnt; i++) xwork[drows[i]] = diag[i] * bwork[drows[i]];
    414:   PetscCall(PetscLogFlops(dcnt));
    415:   PetscCall(VecRestoreArray(red->work, &xwork));
    416:   PetscCall(VecRestoreArrayRead(b, &bwork));
    417:   /* update the right-hand side for the reduced system with diagonal rows (and corresponding columns) removed */
    418:   PetscCall(MatMultTranspose(pc->pmat, x, red->work));
    419:   PetscCall(VecAYPX(red->work, -1.0, b)); /* red->work = b - A^T x */

    421:   PetscCall(VecScatterBegin(red->scatter, red->work, red->b, INSERT_VALUES, SCATTER_FORWARD));
    422:   PetscCall(VecScatterEnd(red->scatter, red->work, red->b, INSERT_VALUES, SCATTER_FORWARD));
    423:   PetscCall(KSPSolveTranspose(red->ksp, red->b, red->x));
    424:   PetscCall(KSPCheckSolve(red->ksp, pc, red->x));
    425:   PetscCall(VecScatterBegin(red->scatter, red->x, x, INSERT_VALUES, SCATTER_REVERSE));
    426:   PetscCall(VecScatterEnd(red->scatter, red->x, x, INSERT_VALUES, SCATTER_REVERSE));
    427:   PetscFunctionReturn(PETSC_SUCCESS);
    428: }

    430: static PetscErrorCode PCDestroy_Redistribute(PC pc)
    431: {
    432:   PC_Redistribute  *red  = (PC_Redistribute *)pc->data;
    433:   PC_FieldSplitLink next = red->splitlinks;

    435:   PetscFunctionBegin;
    436:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetIS_C", NULL));

    438:   while (next) {
    439:     PC_FieldSplitLink ilink;
    440:     PetscCall(PetscFree(next->splitname));
    441:     PetscCall(ISDestroy(&next->is));
    442:     ilink = next;
    443:     next  = next->next;
    444:     PetscCall(PetscFree(ilink));
    445:   }
    446:   PetscCall(VecScatterDestroy(&red->scatter));
    447:   PetscCall(ISDestroy(&red->is));
    448:   PetscCall(VecDestroy(&red->b));
    449:   PetscCall(VecDestroy(&red->x));
    450:   PetscCall(KSPDestroy(&red->ksp));
    451:   PetscCall(VecDestroy(&red->work));
    452:   PetscCall(PetscFree(red->drows));
    453:   PetscCall(PetscFree(red->diag));
    454:   PetscCall(PetscFree(pc->data));
    455:   PetscFunctionReturn(PETSC_SUCCESS);
    456: }

    458: static PetscErrorCode PCSetFromOptions_Redistribute(PC pc, PetscOptionItems *PetscOptionsObject)
    459: {
    460:   PC_Redistribute *red = (PC_Redistribute *)pc->data;

    462:   PetscFunctionBegin;
    463:   PetscCall(KSPSetFromOptions(red->ksp));
    464:   PetscFunctionReturn(PETSC_SUCCESS);
    465: }

    467: /*@
    468:   PCRedistributeGetKSP - Gets the `KSP` created by the `PCREDISTRIBUTE`

    470:   Not Collective

    472:   Input Parameter:
    473: . pc - the preconditioner context

    475:   Output Parameter:
    476: . innerksp - the inner `KSP`

    478:   Level: advanced

    480: .seealso: [](ch_ksp), `KSP`, `PCREDISTRIBUTE`
    481: @*/
    482: PetscErrorCode PCRedistributeGetKSP(PC pc, KSP *innerksp)
    483: {
    484:   PC_Redistribute *red = (PC_Redistribute *)pc->data;

    486:   PetscFunctionBegin;
    488:   PetscAssertPointer(innerksp, 2);
    489:   *innerksp = red->ksp;
    490:   PetscFunctionReturn(PETSC_SUCCESS);
    491: }

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

    497:      Level: intermediate

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

    502:      Usually run this with `-ksp_type preonly`

    504:      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
    505:      -pc_type redistribute -redistribute_ksp_type cg -redistribute_pc_type bjacobi -redistribute_sub_pc_type icc` to take advantage of the symmetry.

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

    510:      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
    511:      between MPI processes inside the preconditioner to balance the number of rows on each process.

    513:      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
    514:      system may be very different (worse) than running that preconditioner on the full system. This is specifically true for elasticity problems.

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

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

    522: PETSC_EXTERN PetscErrorCode PCCreate_Redistribute(PC pc)
    523: {
    524:   PC_Redistribute *red;
    525:   const char      *prefix;

    527:   PetscFunctionBegin;
    528:   PetscCall(PetscNew(&red));
    529:   pc->data = (void *)red;

    531:   pc->ops->apply          = PCApply_Redistribute;
    532:   pc->ops->applytranspose = PCApplyTranspose_Redistribute;
    533:   pc->ops->setup          = PCSetUp_Redistribute;
    534:   pc->ops->destroy        = PCDestroy_Redistribute;
    535:   pc->ops->setfromoptions = PCSetFromOptions_Redistribute;
    536:   pc->ops->view           = PCView_Redistribute;

    538:   PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &red->ksp));
    539:   PetscCall(KSPSetNestLevel(red->ksp, pc->kspnestlevel));
    540:   PetscCall(KSPSetErrorIfNotConverged(red->ksp, pc->erroriffailure));
    541:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)red->ksp, (PetscObject)pc, 1));
    542:   PetscCall(PCGetOptionsPrefix(pc, &prefix));
    543:   PetscCall(KSPSetOptionsPrefix(red->ksp, prefix));
    544:   PetscCall(KSPAppendOptionsPrefix(red->ksp, "redistribute_"));
    545:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetIS_C", PCFieldSplitSetIS_Redistribute));
    546:   PetscFunctionReturn(PETSC_SUCCESS);
    547: }