xref: /petsc/src/ksp/pc/impls/redundant/redundant.c (revision 268865b1601a1129f78d44178bdab76fc555685d)
14b9ad928SBarry Smith /*
23f457be1SHong Zhang   This file defines a "solve the problem redundantly on each subgroup of processor" preconditioner.
34b9ad928SBarry Smith */
4af0996ceSBarry Smith #include <petsc/private/pcimpl.h>
5c6db04a5SJed Brown #include <petscksp.h> /*I "petscksp.h" I*/
64b9ad928SBarry Smith 
74b9ad928SBarry Smith typedef struct {
83e065800SHong Zhang   KSP                ksp;
94b9ad928SBarry Smith   PC                 pc;                    /* actual preconditioner used on each processor */
10ce94432eSBarry Smith   Vec                xsub, ysub;            /* vectors of a subcommunicator to hold parallel vectors of PetscObjectComm((PetscObject)pc) */
113f457be1SHong Zhang   Vec                xdup, ydup;            /* parallel vector that congregates xsub or ysub facilitating vector scattering */
12b3804887SHong Zhang   Mat                pmats;                 /* matrix and optional preconditioner matrix belong to a subcommunicator */
133f457be1SHong Zhang   VecScatter         scatterin, scatterout; /* scatter used to move all values to each processor group (subcommunicator) */
14ace3abfcSBarry Smith   PetscBool          useparallelmat;
15c540e29cSHong Zhang   PetscSubcomm       psubcomm;
161fbd8f88SHong Zhang   PetscInt           nsubcomm; /* num of data structure PetscSubcomm */
17753b7fb9SBarry Smith   PetscBool          shifttypeset;
18753b7fb9SBarry Smith   MatFactorShiftType shifttype;
194b9ad928SBarry Smith } PC_Redundant;
204b9ad928SBarry Smith 
21d71ae5a4SJacob Faibussowitsch PetscErrorCode PCFactorSetShiftType_Redundant(PC pc, MatFactorShiftType shifttype)
22d71ae5a4SJacob Faibussowitsch {
23753b7fb9SBarry Smith   PC_Redundant *red = (PC_Redundant *)pc->data;
24753b7fb9SBarry Smith 
25753b7fb9SBarry Smith   PetscFunctionBegin;
26753b7fb9SBarry Smith   if (red->ksp) {
27753b7fb9SBarry Smith     PC pc;
289566063dSJacob Faibussowitsch     PetscCall(KSPGetPC(red->ksp, &pc));
299566063dSJacob Faibussowitsch     PetscCall(PCFactorSetShiftType(pc, shifttype));
30753b7fb9SBarry Smith   } else {
31753b7fb9SBarry Smith     red->shifttypeset = PETSC_TRUE;
32753b7fb9SBarry Smith     red->shifttype    = shifttype;
33753b7fb9SBarry Smith   }
343ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
35753b7fb9SBarry Smith }
36753b7fb9SBarry Smith 
37d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCView_Redundant(PC pc, PetscViewer viewer)
38d71ae5a4SJacob Faibussowitsch {
394b9ad928SBarry Smith   PC_Redundant *red = (PC_Redundant *)pc->data;
40ace3abfcSBarry Smith   PetscBool     iascii, isstring;
4103ccd0b4SBarry Smith   PetscViewer   subviewer;
424b9ad928SBarry Smith 
434b9ad928SBarry Smith   PetscFunctionBegin;
449566063dSJacob Faibussowitsch   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
459566063dSJacob Faibussowitsch   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
4632077d6dSBarry Smith   if (iascii) {
4703ccd0b4SBarry Smith     if (!red->psubcomm) {
489566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(viewer, "  Not yet setup\n"));
4903ccd0b4SBarry Smith     } else {
5063a3b9bcSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(viewer, "  First (color=0) of %" PetscInt_FMT " PCs follows\n", red->nsubcomm));
519566063dSJacob Faibussowitsch       PetscCall(PetscViewerGetSubViewer(viewer, ((PetscObject)red->pc)->comm, &subviewer));
52f5dd71faSBarry Smith       if (!red->psubcomm->color) { /* only view first redundant pc */
539566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPushTab(subviewer));
549566063dSJacob Faibussowitsch         PetscCall(KSPView(red->ksp, subviewer));
559566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPopTab(subviewer));
564b9ad928SBarry Smith       }
579566063dSJacob Faibussowitsch       PetscCall(PetscViewerRestoreSubViewer(viewer, ((PetscObject)red->pc)->comm, &subviewer));
584b9ad928SBarry Smith     }
5903ccd0b4SBarry Smith   } else if (isstring) {
609566063dSJacob Faibussowitsch     PetscCall(PetscViewerStringSPrintf(viewer, " Redundant solver preconditioner"));
614b9ad928SBarry Smith   }
623ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
634b9ad928SBarry Smith }
644b9ad928SBarry Smith 
6519b3b6edSHong Zhang #include <../src/mat/impls/aij/mpi/mpiaij.h>
66d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCSetUp_Redundant(PC pc)
67d71ae5a4SJacob Faibussowitsch {
684b9ad928SBarry Smith   PC_Redundant *red = (PC_Redundant *)pc->data;
691b81debcSHong Zhang   PetscInt      mstart, mend, mlocal, M;
7013f74950SBarry Smith   PetscMPIInt   size;
71ce94432eSBarry Smith   MPI_Comm      comm, subcomm;
72ddc54837SHong Zhang   Vec           x;
733f457be1SHong Zhang 
744b9ad928SBarry Smith   PetscFunctionBegin;
759566063dSJacob Faibussowitsch   PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
76ddc54837SHong Zhang 
77ddc54837SHong Zhang   /* if pmatrix set by user is sequential then we do not need to gather the parallel matrix */
789566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_size(comm, &size));
79ddc54837SHong Zhang   if (size == 1) red->useparallelmat = PETSC_FALSE;
801fbd8f88SHong Zhang 
814b9ad928SBarry Smith   if (!pc->setupcalled) {
821b81debcSHong Zhang     PetscInt mloc_sub;
8375024027SHong Zhang     if (!red->psubcomm) { /* create red->psubcomm, new ksp and pc over subcomm */
8475024027SHong Zhang       KSP ksp;
859566063dSJacob Faibussowitsch       PetscCall(PCRedundantGetKSP(pc, &ksp));
861b81debcSHong Zhang     }
8775024027SHong Zhang     subcomm = PetscSubcommChild(red->psubcomm);
881fbd8f88SHong Zhang 
891b81debcSHong Zhang     if (red->useparallelmat) {
90*268865b1SPierre Jolivet       /* grab the parallel matrix and put it into the processes of a subcommunicator */
919566063dSJacob Faibussowitsch       PetscCall(MatCreateRedundantMatrix(pc->pmat, red->psubcomm->n, subcomm, MAT_INITIAL_MATRIX, &red->pmats));
92b85f2e9bSHong Zhang 
939566063dSJacob Faibussowitsch       PetscCallMPI(MPI_Comm_size(subcomm, &size));
94b85f2e9bSHong Zhang       if (size > 1) {
9508cecb0aSPierre Jolivet         PetscBool foundpack, issbaij;
969566063dSJacob Faibussowitsch         PetscCall(PetscObjectTypeCompare((PetscObject)red->pmats, MATMPISBAIJ, &issbaij));
9708cecb0aSPierre Jolivet         if (!issbaij) {
989566063dSJacob Faibussowitsch           PetscCall(MatGetFactorAvailable(red->pmats, NULL, MAT_FACTOR_LU, &foundpack));
9908cecb0aSPierre Jolivet         } else {
1009566063dSJacob Faibussowitsch           PetscCall(MatGetFactorAvailable(red->pmats, NULL, MAT_FACTOR_CHOLESKY, &foundpack));
10108cecb0aSPierre Jolivet         }
102b85f2e9bSHong Zhang         if (!foundpack) { /* reset default ksp and pc */
1039566063dSJacob Faibussowitsch           PetscCall(KSPSetType(red->ksp, KSPGMRES));
1049566063dSJacob Faibussowitsch           PetscCall(PCSetType(red->pc, PCBJACOBI));
105c1619fb6SBarry Smith         } else {
1069566063dSJacob Faibussowitsch           PetscCall(PCFactorSetMatSolverType(red->pc, NULL));
107b85f2e9bSHong Zhang         }
108b85f2e9bSHong Zhang       }
109b85f2e9bSHong Zhang 
1109566063dSJacob Faibussowitsch       PetscCall(KSPSetOperators(red->ksp, red->pmats, red->pmats));
1114b9ad928SBarry Smith 
1121b81debcSHong Zhang       /* get working vectors xsub and ysub */
1139566063dSJacob Faibussowitsch       PetscCall(MatCreateVecs(red->pmats, &red->xsub, &red->ysub));
1142fa5cd67SKarl Rupp 
1158b96b0d2SHong Zhang       /* create working vectors xdup and ydup.
1168b96b0d2SHong Zhang        xdup concatenates all xsub's contigously to form a mpi vector over dupcomm  (see PetscSubcommCreate_interlaced())
1178b96b0d2SHong Zhang        ydup concatenates all ysub and has empty local arrays because ysub's arrays will be place into it.
118ce94432eSBarry Smith        Note: we use communicator dupcomm, not PetscObjectComm((PetscObject)pc)! */
1199566063dSJacob Faibussowitsch       PetscCall(MatGetLocalSize(red->pmats, &mloc_sub, NULL));
1209566063dSJacob Faibussowitsch       PetscCall(VecCreateMPI(PetscSubcommContiguousParent(red->psubcomm), mloc_sub, PETSC_DECIDE, &red->xdup));
1219566063dSJacob Faibussowitsch       PetscCall(VecCreateMPIWithArray(PetscSubcommContiguousParent(red->psubcomm), 1, mloc_sub, PETSC_DECIDE, NULL, &red->ydup));
1223f457be1SHong Zhang 
123f68be91cSHong Zhang       /* create vecscatters */
124f68be91cSHong Zhang       if (!red->scatterin) { /* efficiency of scatterin is independent from psubcomm_type! */
1253f457be1SHong Zhang         IS        is1, is2;
1263f457be1SHong Zhang         PetscInt *idx1, *idx2, i, j, k;
12745fc02eaSBarry Smith 
1289566063dSJacob Faibussowitsch         PetscCall(MatCreateVecs(pc->pmat, &x, NULL));
1299566063dSJacob Faibussowitsch         PetscCall(VecGetSize(x, &M));
1309566063dSJacob Faibussowitsch         PetscCall(VecGetOwnershipRange(x, &mstart, &mend));
1311b81debcSHong Zhang         mlocal = mend - mstart;
1329566063dSJacob Faibussowitsch         PetscCall(PetscMalloc2(red->psubcomm->n * mlocal, &idx1, red->psubcomm->n * mlocal, &idx2));
1333f457be1SHong Zhang         j = 0;
1341fbd8f88SHong Zhang         for (k = 0; k < red->psubcomm->n; k++) {
1353f457be1SHong Zhang           for (i = mstart; i < mend; i++) {
1363f457be1SHong Zhang             idx1[j]   = i;
137ddc54837SHong Zhang             idx2[j++] = i + M * k;
1383f457be1SHong Zhang           }
1393f457be1SHong Zhang         }
1409566063dSJacob Faibussowitsch         PetscCall(ISCreateGeneral(comm, red->psubcomm->n * mlocal, idx1, PETSC_COPY_VALUES, &is1));
1419566063dSJacob Faibussowitsch         PetscCall(ISCreateGeneral(comm, red->psubcomm->n * mlocal, idx2, PETSC_COPY_VALUES, &is2));
1429566063dSJacob Faibussowitsch         PetscCall(VecScatterCreate(x, is1, red->xdup, is2, &red->scatterin));
1439566063dSJacob Faibussowitsch         PetscCall(ISDestroy(&is1));
1449566063dSJacob Faibussowitsch         PetscCall(ISDestroy(&is2));
1453f457be1SHong Zhang 
1466909748bSHong Zhang         /* Impl below is good for PETSC_SUBCOMM_INTERLACED (no inter-process communication) and PETSC_SUBCOMM_CONTIGUOUS (communication within subcomm) */
1479566063dSJacob Faibussowitsch         PetscCall(ISCreateStride(comm, mlocal, mstart + red->psubcomm->color * M, 1, &is1));
1489566063dSJacob Faibussowitsch         PetscCall(ISCreateStride(comm, mlocal, mstart, 1, &is2));
1499566063dSJacob Faibussowitsch         PetscCall(VecScatterCreate(red->xdup, is1, x, is2, &red->scatterout));
1509566063dSJacob Faibussowitsch         PetscCall(ISDestroy(&is1));
1519566063dSJacob Faibussowitsch         PetscCall(ISDestroy(&is2));
1529566063dSJacob Faibussowitsch         PetscCall(PetscFree2(idx1, idx2));
1539566063dSJacob Faibussowitsch         PetscCall(VecDestroy(&x));
1541b81debcSHong Zhang       }
155ab661555SHong Zhang     } else { /* !red->useparallelmat */
1569566063dSJacob Faibussowitsch       PetscCall(KSPSetOperators(red->ksp, pc->mat, pc->pmat));
1571b81debcSHong Zhang     }
158ab661555SHong Zhang   } else { /* pc->setupcalled */
1594b9ad928SBarry Smith     if (red->useparallelmat) {
160ab661555SHong Zhang       MatReuse reuse;
161*268865b1SPierre Jolivet       /* grab the parallel matrix and put it into the processes of a subcommunicator */
1621b81debcSHong Zhang       /*--------------------------------------------------------------------------*/
163ab661555SHong Zhang       if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1644b9ad928SBarry Smith         /* destroy old matrices */
1659566063dSJacob Faibussowitsch         PetscCall(MatDestroy(&red->pmats));
166ab661555SHong Zhang         reuse = MAT_INITIAL_MATRIX;
1674b9ad928SBarry Smith       } else {
168ab661555SHong Zhang         reuse = MAT_REUSE_MATRIX;
169ab661555SHong Zhang       }
1709566063dSJacob Faibussowitsch       PetscCall(MatCreateRedundantMatrix(pc->pmat, red->psubcomm->n, PetscSubcommChild(red->psubcomm), reuse, &red->pmats));
1719566063dSJacob Faibussowitsch       PetscCall(KSPSetOperators(red->ksp, red->pmats, red->pmats));
172ab661555SHong Zhang     } else { /* !red->useparallelmat */
1739566063dSJacob Faibussowitsch       PetscCall(KSPSetOperators(red->ksp, pc->mat, pc->pmat));
1744b9ad928SBarry Smith     }
175ab661555SHong Zhang   }
1761b81debcSHong Zhang 
1771baa6e33SBarry Smith   if (pc->setfromoptionscalled) PetscCall(KSPSetFromOptions(red->ksp));
1789566063dSJacob Faibussowitsch   PetscCall(KSPSetUp(red->ksp));
1793ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
1804b9ad928SBarry Smith }
1814b9ad928SBarry Smith 
182d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCApply_Redundant(PC pc, Vec x, Vec y)
183d71ae5a4SJacob Faibussowitsch {
1844b9ad928SBarry Smith   PC_Redundant *red = (PC_Redundant *)pc->data;
1853f457be1SHong Zhang   PetscScalar  *array;
1864b9ad928SBarry Smith 
1874b9ad928SBarry Smith   PetscFunctionBegin;
188ddc54837SHong Zhang   if (!red->useparallelmat) {
1899566063dSJacob Faibussowitsch     PetscCall(KSPSolve(red->ksp, x, y));
1909566063dSJacob Faibussowitsch     PetscCall(KSPCheckSolve(red->ksp, pc, y));
1913ba16761SJacob Faibussowitsch     PetscFunctionReturn(PETSC_SUCCESS);
192ddc54837SHong Zhang   }
193ddc54837SHong Zhang 
1943f457be1SHong Zhang   /* scatter x to xdup */
1959566063dSJacob Faibussowitsch   PetscCall(VecScatterBegin(red->scatterin, x, red->xdup, INSERT_VALUES, SCATTER_FORWARD));
1969566063dSJacob Faibussowitsch   PetscCall(VecScatterEnd(red->scatterin, x, red->xdup, INSERT_VALUES, SCATTER_FORWARD));
1973f457be1SHong Zhang 
1983f457be1SHong Zhang   /* place xdup's local array into xsub */
1999566063dSJacob Faibussowitsch   PetscCall(VecGetArray(red->xdup, &array));
2009566063dSJacob Faibussowitsch   PetscCall(VecPlaceArray(red->xsub, (const PetscScalar *)array));
2014b9ad928SBarry Smith 
2024b9ad928SBarry Smith   /* apply preconditioner on each processor */
2039566063dSJacob Faibussowitsch   PetscCall(KSPSolve(red->ksp, red->xsub, red->ysub));
2049566063dSJacob Faibussowitsch   PetscCall(KSPCheckSolve(red->ksp, pc, red->ysub));
2059566063dSJacob Faibussowitsch   PetscCall(VecResetArray(red->xsub));
2069566063dSJacob Faibussowitsch   PetscCall(VecRestoreArray(red->xdup, &array));
2074b9ad928SBarry Smith 
2083f457be1SHong Zhang   /* place ysub's local array into ydup */
2099566063dSJacob Faibussowitsch   PetscCall(VecGetArray(red->ysub, &array));
2109566063dSJacob Faibussowitsch   PetscCall(VecPlaceArray(red->ydup, (const PetscScalar *)array));
2113f457be1SHong Zhang 
2123f457be1SHong Zhang   /* scatter ydup to y */
2139566063dSJacob Faibussowitsch   PetscCall(VecScatterBegin(red->scatterout, red->ydup, y, INSERT_VALUES, SCATTER_FORWARD));
2149566063dSJacob Faibussowitsch   PetscCall(VecScatterEnd(red->scatterout, red->ydup, y, INSERT_VALUES, SCATTER_FORWARD));
2159566063dSJacob Faibussowitsch   PetscCall(VecResetArray(red->ydup));
2169566063dSJacob Faibussowitsch   PetscCall(VecRestoreArray(red->ysub, &array));
2173ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
2184b9ad928SBarry Smith }
2194b9ad928SBarry Smith 
220d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCApplyTranspose_Redundant(PC pc, Vec x, Vec y)
221d71ae5a4SJacob Faibussowitsch {
222d88bfacbSStefano Zampini   PC_Redundant *red = (PC_Redundant *)pc->data;
223d88bfacbSStefano Zampini   PetscScalar  *array;
224d88bfacbSStefano Zampini 
225d88bfacbSStefano Zampini   PetscFunctionBegin;
226d88bfacbSStefano Zampini   if (!red->useparallelmat) {
2279566063dSJacob Faibussowitsch     PetscCall(KSPSolveTranspose(red->ksp, x, y));
2289566063dSJacob Faibussowitsch     PetscCall(KSPCheckSolve(red->ksp, pc, y));
2293ba16761SJacob Faibussowitsch     PetscFunctionReturn(PETSC_SUCCESS);
230d88bfacbSStefano Zampini   }
231d88bfacbSStefano Zampini 
232d88bfacbSStefano Zampini   /* scatter x to xdup */
2339566063dSJacob Faibussowitsch   PetscCall(VecScatterBegin(red->scatterin, x, red->xdup, INSERT_VALUES, SCATTER_FORWARD));
2349566063dSJacob Faibussowitsch   PetscCall(VecScatterEnd(red->scatterin, x, red->xdup, INSERT_VALUES, SCATTER_FORWARD));
235d88bfacbSStefano Zampini 
236d88bfacbSStefano Zampini   /* place xdup's local array into xsub */
2379566063dSJacob Faibussowitsch   PetscCall(VecGetArray(red->xdup, &array));
2389566063dSJacob Faibussowitsch   PetscCall(VecPlaceArray(red->xsub, (const PetscScalar *)array));
239d88bfacbSStefano Zampini 
240d88bfacbSStefano Zampini   /* apply preconditioner on each processor */
2419566063dSJacob Faibussowitsch   PetscCall(KSPSolveTranspose(red->ksp, red->xsub, red->ysub));
2429566063dSJacob Faibussowitsch   PetscCall(KSPCheckSolve(red->ksp, pc, red->ysub));
2439566063dSJacob Faibussowitsch   PetscCall(VecResetArray(red->xsub));
2449566063dSJacob Faibussowitsch   PetscCall(VecRestoreArray(red->xdup, &array));
245d88bfacbSStefano Zampini 
246d88bfacbSStefano Zampini   /* place ysub's local array into ydup */
2479566063dSJacob Faibussowitsch   PetscCall(VecGetArray(red->ysub, &array));
2489566063dSJacob Faibussowitsch   PetscCall(VecPlaceArray(red->ydup, (const PetscScalar *)array));
249d88bfacbSStefano Zampini 
250d88bfacbSStefano Zampini   /* scatter ydup to y */
2519566063dSJacob Faibussowitsch   PetscCall(VecScatterBegin(red->scatterout, red->ydup, y, INSERT_VALUES, SCATTER_FORWARD));
2529566063dSJacob Faibussowitsch   PetscCall(VecScatterEnd(red->scatterout, red->ydup, y, INSERT_VALUES, SCATTER_FORWARD));
2539566063dSJacob Faibussowitsch   PetscCall(VecResetArray(red->ydup));
2549566063dSJacob Faibussowitsch   PetscCall(VecRestoreArray(red->ysub, &array));
2553ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
256d88bfacbSStefano Zampini }
257d88bfacbSStefano Zampini 
258d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCReset_Redundant(PC pc)
259d71ae5a4SJacob Faibussowitsch {
2604b9ad928SBarry Smith   PC_Redundant *red = (PC_Redundant *)pc->data;
2614b9ad928SBarry Smith 
2624b9ad928SBarry Smith   PetscFunctionBegin;
2631b81debcSHong Zhang   if (red->useparallelmat) {
2649566063dSJacob Faibussowitsch     PetscCall(VecScatterDestroy(&red->scatterin));
2659566063dSJacob Faibussowitsch     PetscCall(VecScatterDestroy(&red->scatterout));
2669566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&red->ysub));
2679566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&red->xsub));
2689566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&red->xdup));
2699566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&red->ydup));
2701b81debcSHong Zhang   }
2719566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&red->pmats));
2729566063dSJacob Faibussowitsch   PetscCall(KSPReset(red->ksp));
2733ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
2741ea5a559SBarry Smith }
2751ea5a559SBarry Smith 
276d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCDestroy_Redundant(PC pc)
277d71ae5a4SJacob Faibussowitsch {
2781ea5a559SBarry Smith   PC_Redundant *red = (PC_Redundant *)pc->data;
2791ea5a559SBarry Smith 
2801ea5a559SBarry Smith   PetscFunctionBegin;
2819566063dSJacob Faibussowitsch   PetscCall(PCReset_Redundant(pc));
2829566063dSJacob Faibussowitsch   PetscCall(KSPDestroy(&red->ksp));
2839566063dSJacob Faibussowitsch   PetscCall(PetscSubcommDestroy(&red->psubcomm));
2842e956fe4SStefano Zampini   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantSetScatter_C", NULL));
2852e956fe4SStefano Zampini   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantSetNumber_C", NULL));
2862e956fe4SStefano Zampini   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantGetKSP_C", NULL));
2872e956fe4SStefano Zampini   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantGetOperators_C", NULL));
2882e956fe4SStefano Zampini   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFactorSetShiftType_C", NULL));
2899566063dSJacob Faibussowitsch   PetscCall(PetscFree(pc->data));
2903ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
2914b9ad928SBarry Smith }
2924b9ad928SBarry Smith 
293d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCSetFromOptions_Redundant(PC pc, PetscOptionItems *PetscOptionsObject)
294d71ae5a4SJacob Faibussowitsch {
295a98ce0f4SHong Zhang   PC_Redundant *red = (PC_Redundant *)pc->data;
296a98ce0f4SHong Zhang 
2974b9ad928SBarry Smith   PetscFunctionBegin;
298d0609cedSBarry Smith   PetscOptionsHeadBegin(PetscOptionsObject, "Redundant options");
2999566063dSJacob Faibussowitsch   PetscCall(PetscOptionsInt("-pc_redundant_number", "Number of redundant pc", "PCRedundantSetNumber", red->nsubcomm, &red->nsubcomm, NULL));
300d0609cedSBarry Smith   PetscOptionsHeadEnd();
3013ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
3024b9ad928SBarry Smith }
3034b9ad928SBarry Smith 
304d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCRedundantSetNumber_Redundant(PC pc, PetscInt nreds)
305d71ae5a4SJacob Faibussowitsch {
30609a6bc64SHong Zhang   PC_Redundant *red = (PC_Redundant *)pc->data;
30709a6bc64SHong Zhang 
30809a6bc64SHong Zhang   PetscFunctionBegin;
30909a6bc64SHong Zhang   red->nsubcomm = nreds;
3103ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
31109a6bc64SHong Zhang }
31209a6bc64SHong Zhang 
31309a6bc64SHong Zhang /*@
31409a6bc64SHong Zhang    PCRedundantSetNumber - Sets the number of redundant preconditioner contexts.
31509a6bc64SHong Zhang 
316c3339decSBarry Smith    Logically Collective
31709a6bc64SHong Zhang 
31809a6bc64SHong Zhang    Input Parameters:
31909a6bc64SHong Zhang +  pc - the preconditioner context
3209b21d695SBarry Smith -  nredundant - number of redundant preconditioner contexts; for example if you are using 64 MPI processes and
3219b21d695SBarry Smith                               use an nredundant of 4 there will be 4 parallel solves each on 16 = 64/4 processes.
32209a6bc64SHong Zhang 
32309a6bc64SHong Zhang    Level: advanced
32409a6bc64SHong Zhang 
325f1580f4eSBarry Smith .seealso: `PCREDUNDANT`
32609a6bc64SHong Zhang @*/
327d71ae5a4SJacob Faibussowitsch PetscErrorCode PCRedundantSetNumber(PC pc, PetscInt nredundant)
328d71ae5a4SJacob Faibussowitsch {
32909a6bc64SHong Zhang   PetscFunctionBegin;
3300700a824SBarry Smith   PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
33163a3b9bcSJacob Faibussowitsch   PetscCheck(nredundant > 0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "num of redundant pc %" PetscInt_FMT " must be positive", nredundant);
332cac4c232SBarry Smith   PetscTryMethod(pc, "PCRedundantSetNumber_C", (PC, PetscInt), (pc, nredundant));
3333ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
33409a6bc64SHong Zhang }
33509a6bc64SHong Zhang 
336d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCRedundantSetScatter_Redundant(PC pc, VecScatter in, VecScatter out)
337d71ae5a4SJacob Faibussowitsch {
3384b9ad928SBarry Smith   PC_Redundant *red = (PC_Redundant *)pc->data;
3394b9ad928SBarry Smith 
3404b9ad928SBarry Smith   PetscFunctionBegin;
3419566063dSJacob Faibussowitsch   PetscCall(PetscObjectReference((PetscObject)in));
3429566063dSJacob Faibussowitsch   PetscCall(VecScatterDestroy(&red->scatterin));
3432fa5cd67SKarl Rupp 
344c3122656SLisandro Dalcin   red->scatterin = in;
3452fa5cd67SKarl Rupp 
3469566063dSJacob Faibussowitsch   PetscCall(PetscObjectReference((PetscObject)out));
3479566063dSJacob Faibussowitsch   PetscCall(VecScatterDestroy(&red->scatterout));
348c3122656SLisandro Dalcin   red->scatterout = out;
3493ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
3504b9ad928SBarry Smith }
3514b9ad928SBarry Smith 
3524b9ad928SBarry Smith /*@
3534b9ad928SBarry Smith    PCRedundantSetScatter - Sets the scatter used to copy values into the
3544b9ad928SBarry Smith      redundant local solve and the scatter to move them back into the global
3554b9ad928SBarry Smith      vector.
3564b9ad928SBarry Smith 
357c3339decSBarry Smith    Logically Collective
3584b9ad928SBarry Smith 
3594b9ad928SBarry Smith    Input Parameters:
3604b9ad928SBarry Smith +  pc - the preconditioner context
3614b9ad928SBarry Smith .  in - the scatter to move the values in
3624b9ad928SBarry Smith -  out - the scatter to move them out
3634b9ad928SBarry Smith 
3644b9ad928SBarry Smith    Level: advanced
3654b9ad928SBarry Smith 
366f1580f4eSBarry Smith .seealso: `PCREDUNDANT`
3674b9ad928SBarry Smith @*/
368d71ae5a4SJacob Faibussowitsch PetscErrorCode PCRedundantSetScatter(PC pc, VecScatter in, VecScatter out)
369d71ae5a4SJacob Faibussowitsch {
3704b9ad928SBarry Smith   PetscFunctionBegin;
3710700a824SBarry Smith   PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
37297929ea7SJunchao Zhang   PetscValidHeaderSpecific(in, PETSCSF_CLASSID, 2);
37397929ea7SJunchao Zhang   PetscValidHeaderSpecific(out, PETSCSF_CLASSID, 3);
374cac4c232SBarry Smith   PetscTryMethod(pc, "PCRedundantSetScatter_C", (PC, VecScatter, VecScatter), (pc, in, out));
3753ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
3764b9ad928SBarry Smith }
3774b9ad928SBarry Smith 
378d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCRedundantGetKSP_Redundant(PC pc, KSP *innerksp)
379d71ae5a4SJacob Faibussowitsch {
3804b9ad928SBarry Smith   PC_Redundant *red = (PC_Redundant *)pc->data;
38175024027SHong Zhang   MPI_Comm      comm, subcomm;
38275024027SHong Zhang   const char   *prefix;
38308cecb0aSPierre Jolivet   PetscBool     issbaij;
3844b9ad928SBarry Smith 
3854b9ad928SBarry Smith   PetscFunctionBegin;
38675024027SHong Zhang   if (!red->psubcomm) {
3879566063dSJacob Faibussowitsch     PetscCall(PCGetOptionsPrefix(pc, &prefix));
388e5acf8a4SHong Zhang 
3899566063dSJacob Faibussowitsch     PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
3909566063dSJacob Faibussowitsch     PetscCall(PetscSubcommCreate(comm, &red->psubcomm));
3919566063dSJacob Faibussowitsch     PetscCall(PetscSubcommSetNumber(red->psubcomm, red->nsubcomm));
3929566063dSJacob Faibussowitsch     PetscCall(PetscSubcommSetType(red->psubcomm, PETSC_SUBCOMM_CONTIGUOUS));
393e5acf8a4SHong Zhang 
3949566063dSJacob Faibussowitsch     PetscCall(PetscSubcommSetOptionsPrefix(red->psubcomm, prefix));
3959566063dSJacob Faibussowitsch     PetscCall(PetscSubcommSetFromOptions(red->psubcomm));
39675024027SHong Zhang 
39775024027SHong Zhang     /* create a new PC that processors in each subcomm have copy of */
39875024027SHong Zhang     subcomm = PetscSubcommChild(red->psubcomm);
39975024027SHong Zhang 
4009566063dSJacob Faibussowitsch     PetscCall(KSPCreate(subcomm, &red->ksp));
4019566063dSJacob Faibussowitsch     PetscCall(KSPSetErrorIfNotConverged(red->ksp, pc->erroriffailure));
4029566063dSJacob Faibussowitsch     PetscCall(PetscObjectIncrementTabLevel((PetscObject)red->ksp, (PetscObject)pc, 1));
4039566063dSJacob Faibussowitsch     PetscCall(KSPSetType(red->ksp, KSPPREONLY));
4049566063dSJacob Faibussowitsch     PetscCall(KSPGetPC(red->ksp, &red->pc));
4059566063dSJacob Faibussowitsch     PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATSEQSBAIJ, &issbaij));
40648a46eb9SPierre Jolivet     if (!issbaij) PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATMPISBAIJ, &issbaij));
40708cecb0aSPierre Jolivet     if (!issbaij) {
4089566063dSJacob Faibussowitsch       PetscCall(PCSetType(red->pc, PCLU));
40908cecb0aSPierre Jolivet     } else {
4109566063dSJacob Faibussowitsch       PetscCall(PCSetType(red->pc, PCCHOLESKY));
41108cecb0aSPierre Jolivet     }
412753b7fb9SBarry Smith     if (red->shifttypeset) {
4139566063dSJacob Faibussowitsch       PetscCall(PCFactorSetShiftType(red->pc, red->shifttype));
414753b7fb9SBarry Smith       red->shifttypeset = PETSC_FALSE;
415753b7fb9SBarry Smith     }
4169566063dSJacob Faibussowitsch     PetscCall(KSPSetOptionsPrefix(red->ksp, prefix));
4179566063dSJacob Faibussowitsch     PetscCall(KSPAppendOptionsPrefix(red->ksp, "redundant_"));
41875024027SHong Zhang   }
41983ab6a24SBarry Smith   *innerksp = red->ksp;
4203ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
4214b9ad928SBarry Smith }
4224b9ad928SBarry Smith 
4234b9ad928SBarry Smith /*@
424f1580f4eSBarry Smith    PCRedundantGetKSP - Gets the less parallel `KSP` created by the redundant `PC`.
4254b9ad928SBarry Smith 
4264b9ad928SBarry Smith    Not Collective
4274b9ad928SBarry Smith 
4284b9ad928SBarry Smith    Input Parameter:
4294b9ad928SBarry Smith .  pc - the preconditioner context
4304b9ad928SBarry Smith 
4314b9ad928SBarry Smith    Output Parameter:
432f1580f4eSBarry Smith .  innerksp - the `KSP` on the smaller set of processes
4334b9ad928SBarry Smith 
4344b9ad928SBarry Smith    Level: advanced
4354b9ad928SBarry Smith 
436f1580f4eSBarry Smith .seealso: `PCREDUNDANT`
4374b9ad928SBarry Smith @*/
438d71ae5a4SJacob Faibussowitsch PetscErrorCode PCRedundantGetKSP(PC pc, KSP *innerksp)
439d71ae5a4SJacob Faibussowitsch {
4404b9ad928SBarry Smith   PetscFunctionBegin;
4410700a824SBarry Smith   PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
44283ab6a24SBarry Smith   PetscValidPointer(innerksp, 2);
443cac4c232SBarry Smith   PetscUseMethod(pc, "PCRedundantGetKSP_C", (PC, KSP *), (pc, innerksp));
4443ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
4454b9ad928SBarry Smith }
4464b9ad928SBarry Smith 
447d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCRedundantGetOperators_Redundant(PC pc, Mat *mat, Mat *pmat)
448d71ae5a4SJacob Faibussowitsch {
4494b9ad928SBarry Smith   PC_Redundant *red = (PC_Redundant *)pc->data;
4504b9ad928SBarry Smith 
4514b9ad928SBarry Smith   PetscFunctionBegin;
452b3804887SHong Zhang   if (mat) *mat = red->pmats;
453b3804887SHong Zhang   if (pmat) *pmat = red->pmats;
4543ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
4554b9ad928SBarry Smith }
4564b9ad928SBarry Smith 
4574b9ad928SBarry Smith /*@
4584b9ad928SBarry Smith    PCRedundantGetOperators - gets the sequential matrix and preconditioner matrix
4594b9ad928SBarry Smith 
4604b9ad928SBarry Smith    Not Collective
4614b9ad928SBarry Smith 
4624b9ad928SBarry Smith    Input Parameter:
4634b9ad928SBarry Smith .  pc - the preconditioner context
4644b9ad928SBarry Smith 
4654b9ad928SBarry Smith    Output Parameters:
4664b9ad928SBarry Smith +  mat - the matrix
4674b9ad928SBarry Smith -  pmat - the (possibly different) preconditioner matrix
4684b9ad928SBarry Smith 
4694b9ad928SBarry Smith    Level: advanced
4704b9ad928SBarry Smith 
471f1580f4eSBarry Smith .seealso: `PCREDUNDANT`
4724b9ad928SBarry Smith @*/
473d71ae5a4SJacob Faibussowitsch PetscErrorCode PCRedundantGetOperators(PC pc, Mat *mat, Mat *pmat)
474d71ae5a4SJacob Faibussowitsch {
4754b9ad928SBarry Smith   PetscFunctionBegin;
4760700a824SBarry Smith   PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
4774482741eSBarry Smith   if (mat) PetscValidPointer(mat, 2);
4784482741eSBarry Smith   if (pmat) PetscValidPointer(pmat, 3);
479cac4c232SBarry Smith   PetscUseMethod(pc, "PCRedundantGetOperators_C", (PC, Mat *, Mat *), (pc, mat, pmat));
4803ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
4814b9ad928SBarry Smith }
4824b9ad928SBarry Smith 
48337a17b4dSBarry Smith /*MC
484f1580f4eSBarry Smith      PCREDUNDANT - Runs a `KSP` solver with preconditioner for the entire problem on subgroups of processors
48537a17b4dSBarry Smith 
486f1580f4eSBarry Smith   Options Database Key:
4879b21d695SBarry Smith .  -pc_redundant_number <n> - number of redundant solves, for example if you are using 64 MPI processes and
4889b21d695SBarry Smith                               use an n of 4 there will be 4 parallel solves each on 16 = 64/4 processes.
48909391456SBarry Smith 
49037a17b4dSBarry Smith    Level: intermediate
49137a17b4dSBarry Smith 
49295452b02SPatrick Sanan    Notes:
493f1580f4eSBarry Smith    Options for the redundant preconditioners can be set using the options database prefix -redundant_
49483ab6a24SBarry Smith 
495f1580f4eSBarry Smith    The default `KSP` is preonly and the default `PC` is `PCLU` or `PCCHOLESKY` if Pmat is of type `MATSBAIJ`.
496753b7fb9SBarry Smith 
497f1580f4eSBarry Smith    `PCFactorSetShiftType()` applied to this `PC` will convey they shift type into the inner `PC` if it is factorization based.
498f1580f4eSBarry Smith 
499f1580f4eSBarry Smith    Developer Note:
500f1580f4eSBarry Smith    `PCSetInitialGuessNonzero()` is not used by this class but likely should be.
5019cfaa89bSBarry Smith 
502db781477SPatrick Sanan .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PCRedundantSetScatter()`,
503f1580f4eSBarry Smith           `PCRedundantGetKSP()`, `PCRedundantGetOperators()`, `PCRedundantSetNumber()`, `PCREDISTRIBUTE`
50437a17b4dSBarry Smith M*/
50537a17b4dSBarry Smith 
506d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode PCCreate_Redundant(PC pc)
507d71ae5a4SJacob Faibussowitsch {
5084b9ad928SBarry Smith   PC_Redundant *red;
50969db28dcSHong Zhang   PetscMPIInt   size;
5103f457be1SHong Zhang 
5114b9ad928SBarry Smith   PetscFunctionBegin;
5124dfa11a4SJacob Faibussowitsch   PetscCall(PetscNew(&red));
5139566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
5142fa5cd67SKarl Rupp 
51569db28dcSHong Zhang   red->nsubcomm       = size;
5164b9ad928SBarry Smith   red->useparallelmat = PETSC_TRUE;
5171fbd8f88SHong Zhang   pc->data            = (void *)red;
5184b9ad928SBarry Smith 
5194b9ad928SBarry Smith   pc->ops->apply          = PCApply_Redundant;
520d88bfacbSStefano Zampini   pc->ops->applytranspose = PCApplyTranspose_Redundant;
5214b9ad928SBarry Smith   pc->ops->setup          = PCSetUp_Redundant;
5224b9ad928SBarry Smith   pc->ops->destroy        = PCDestroy_Redundant;
5231ea5a559SBarry Smith   pc->ops->reset          = PCReset_Redundant;
5244b9ad928SBarry Smith   pc->ops->setfromoptions = PCSetFromOptions_Redundant;
5254b9ad928SBarry Smith   pc->ops->view           = PCView_Redundant;
5262fa5cd67SKarl Rupp 
5279566063dSJacob Faibussowitsch   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantSetScatter_C", PCRedundantSetScatter_Redundant));
5289566063dSJacob Faibussowitsch   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantSetNumber_C", PCRedundantSetNumber_Redundant));
5299566063dSJacob Faibussowitsch   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantGetKSP_C", PCRedundantGetKSP_Redundant));
5309566063dSJacob Faibussowitsch   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantGetOperators_C", PCRedundantGetOperators_Redundant));
5319566063dSJacob Faibussowitsch   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFactorSetShiftType_C", PCFactorSetShiftType_Redundant));
5323ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
5334b9ad928SBarry Smith }
534