xref: /petsc/src/ksp/pc/impls/redundant/redundant.c (revision 4dfa11a44d5adf2389f1d3acbc8f3c1116dc6c3a)
1dba47a55SKris Buschelman 
24b9ad928SBarry Smith /*
33f457be1SHong Zhang   This file defines a "solve the problem redundantly on each subgroup of processor" preconditioner.
44b9ad928SBarry Smith */
5af0996ceSBarry Smith #include <petsc/private/pcimpl.h>
6c6db04a5SJed Brown #include <petscksp.h> /*I "petscksp.h" I*/
74b9ad928SBarry Smith 
84b9ad928SBarry Smith typedef struct {
93e065800SHong Zhang   KSP                ksp;
104b9ad928SBarry Smith   PC                 pc;                    /* actual preconditioner used on each processor */
11ce94432eSBarry Smith   Vec                xsub, ysub;            /* vectors of a subcommunicator to hold parallel vectors of PetscObjectComm((PetscObject)pc) */
123f457be1SHong Zhang   Vec                xdup, ydup;            /* parallel vector that congregates xsub or ysub facilitating vector scattering */
13b3804887SHong Zhang   Mat                pmats;                 /* matrix and optional preconditioner matrix belong to a subcommunicator */
143f457be1SHong Zhang   VecScatter         scatterin, scatterout; /* scatter used to move all values to each processor group (subcommunicator) */
15ace3abfcSBarry Smith   PetscBool          useparallelmat;
16c540e29cSHong Zhang   PetscSubcomm       psubcomm;
171fbd8f88SHong Zhang   PetscInt           nsubcomm; /* num of data structure PetscSubcomm */
18753b7fb9SBarry Smith   PetscBool          shifttypeset;
19753b7fb9SBarry Smith   MatFactorShiftType shifttype;
204b9ad928SBarry Smith } PC_Redundant;
214b9ad928SBarry Smith 
229371c9d4SSatish Balay PetscErrorCode PCFactorSetShiftType_Redundant(PC pc, MatFactorShiftType shifttype) {
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   }
34753b7fb9SBarry Smith   PetscFunctionReturn(0);
35753b7fb9SBarry Smith }
36753b7fb9SBarry Smith 
379371c9d4SSatish Balay static PetscErrorCode PCView_Redundant(PC pc, PetscViewer viewer) {
384b9ad928SBarry Smith   PC_Redundant *red = (PC_Redundant *)pc->data;
39ace3abfcSBarry Smith   PetscBool     iascii, isstring;
4003ccd0b4SBarry Smith   PetscViewer   subviewer;
414b9ad928SBarry Smith 
424b9ad928SBarry Smith   PetscFunctionBegin;
439566063dSJacob Faibussowitsch   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
449566063dSJacob Faibussowitsch   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
4532077d6dSBarry Smith   if (iascii) {
4603ccd0b4SBarry Smith     if (!red->psubcomm) {
479566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(viewer, "  Not yet setup\n"));
4803ccd0b4SBarry Smith     } else {
4963a3b9bcSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(viewer, "  First (color=0) of %" PetscInt_FMT " PCs follows\n", red->nsubcomm));
509566063dSJacob Faibussowitsch       PetscCall(PetscViewerGetSubViewer(viewer, ((PetscObject)red->pc)->comm, &subviewer));
51f5dd71faSBarry Smith       if (!red->psubcomm->color) { /* only view first redundant pc */
529566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPushTab(subviewer));
539566063dSJacob Faibussowitsch         PetscCall(KSPView(red->ksp, subviewer));
549566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPopTab(subviewer));
554b9ad928SBarry Smith       }
569566063dSJacob Faibussowitsch       PetscCall(PetscViewerRestoreSubViewer(viewer, ((PetscObject)red->pc)->comm, &subviewer));
574b9ad928SBarry Smith     }
5803ccd0b4SBarry Smith   } else if (isstring) {
599566063dSJacob Faibussowitsch     PetscCall(PetscViewerStringSPrintf(viewer, " Redundant solver preconditioner"));
604b9ad928SBarry Smith   }
614b9ad928SBarry Smith   PetscFunctionReturn(0);
624b9ad928SBarry Smith }
634b9ad928SBarry Smith 
6419b3b6edSHong Zhang #include <../src/mat/impls/aij/mpi/mpiaij.h>
659371c9d4SSatish Balay static PetscErrorCode PCSetUp_Redundant(PC pc) {
664b9ad928SBarry Smith   PC_Redundant *red = (PC_Redundant *)pc->data;
671b81debcSHong Zhang   PetscInt      mstart, mend, mlocal, M;
6813f74950SBarry Smith   PetscMPIInt   size;
69ce94432eSBarry Smith   MPI_Comm      comm, subcomm;
70ddc54837SHong Zhang   Vec           x;
713f457be1SHong Zhang 
724b9ad928SBarry Smith   PetscFunctionBegin;
739566063dSJacob Faibussowitsch   PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
74ddc54837SHong Zhang 
75ddc54837SHong Zhang   /* if pmatrix set by user is sequential then we do not need to gather the parallel matrix */
769566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_size(comm, &size));
77ddc54837SHong Zhang   if (size == 1) red->useparallelmat = PETSC_FALSE;
781fbd8f88SHong Zhang 
794b9ad928SBarry Smith   if (!pc->setupcalled) {
801b81debcSHong Zhang     PetscInt mloc_sub;
8175024027SHong Zhang     if (!red->psubcomm) { /* create red->psubcomm, new ksp and pc over subcomm */
8275024027SHong Zhang       KSP ksp;
839566063dSJacob Faibussowitsch       PetscCall(PCRedundantGetKSP(pc, &ksp));
841b81debcSHong Zhang     }
8575024027SHong Zhang     subcomm = PetscSubcommChild(red->psubcomm);
861fbd8f88SHong Zhang 
871b81debcSHong Zhang     if (red->useparallelmat) {
881b81debcSHong Zhang       /* grab the parallel matrix and put it into processors of a subcomminicator */
899566063dSJacob Faibussowitsch       PetscCall(MatCreateRedundantMatrix(pc->pmat, red->psubcomm->n, subcomm, MAT_INITIAL_MATRIX, &red->pmats));
90b85f2e9bSHong Zhang 
919566063dSJacob Faibussowitsch       PetscCallMPI(MPI_Comm_size(subcomm, &size));
92b85f2e9bSHong Zhang       if (size > 1) {
9308cecb0aSPierre Jolivet         PetscBool foundpack, issbaij;
949566063dSJacob Faibussowitsch         PetscCall(PetscObjectTypeCompare((PetscObject)red->pmats, MATMPISBAIJ, &issbaij));
9508cecb0aSPierre Jolivet         if (!issbaij) {
969566063dSJacob Faibussowitsch           PetscCall(MatGetFactorAvailable(red->pmats, NULL, MAT_FACTOR_LU, &foundpack));
9708cecb0aSPierre Jolivet         } else {
989566063dSJacob Faibussowitsch           PetscCall(MatGetFactorAvailable(red->pmats, NULL, MAT_FACTOR_CHOLESKY, &foundpack));
9908cecb0aSPierre Jolivet         }
100b85f2e9bSHong Zhang         if (!foundpack) { /* reset default ksp and pc */
1019566063dSJacob Faibussowitsch           PetscCall(KSPSetType(red->ksp, KSPGMRES));
1029566063dSJacob Faibussowitsch           PetscCall(PCSetType(red->pc, PCBJACOBI));
103c1619fb6SBarry Smith         } else {
1049566063dSJacob Faibussowitsch           PetscCall(PCFactorSetMatSolverType(red->pc, NULL));
105b85f2e9bSHong Zhang         }
106b85f2e9bSHong Zhang       }
107b85f2e9bSHong Zhang 
1089566063dSJacob Faibussowitsch       PetscCall(KSPSetOperators(red->ksp, red->pmats, red->pmats));
1094b9ad928SBarry Smith 
1101b81debcSHong Zhang       /* get working vectors xsub and ysub */
1119566063dSJacob Faibussowitsch       PetscCall(MatCreateVecs(red->pmats, &red->xsub, &red->ysub));
1122fa5cd67SKarl Rupp 
1138b96b0d2SHong Zhang       /* create working vectors xdup and ydup.
1148b96b0d2SHong Zhang        xdup concatenates all xsub's contigously to form a mpi vector over dupcomm  (see PetscSubcommCreate_interlaced())
1158b96b0d2SHong Zhang        ydup concatenates all ysub and has empty local arrays because ysub's arrays will be place into it.
116ce94432eSBarry Smith        Note: we use communicator dupcomm, not PetscObjectComm((PetscObject)pc)! */
1179566063dSJacob Faibussowitsch       PetscCall(MatGetLocalSize(red->pmats, &mloc_sub, NULL));
1189566063dSJacob Faibussowitsch       PetscCall(VecCreateMPI(PetscSubcommContiguousParent(red->psubcomm), mloc_sub, PETSC_DECIDE, &red->xdup));
1199566063dSJacob Faibussowitsch       PetscCall(VecCreateMPIWithArray(PetscSubcommContiguousParent(red->psubcomm), 1, mloc_sub, PETSC_DECIDE, NULL, &red->ydup));
1203f457be1SHong Zhang 
121f68be91cSHong Zhang       /* create vecscatters */
122f68be91cSHong Zhang       if (!red->scatterin) { /* efficiency of scatterin is independent from psubcomm_type! */
1233f457be1SHong Zhang         IS        is1, is2;
1243f457be1SHong Zhang         PetscInt *idx1, *idx2, i, j, k;
12545fc02eaSBarry Smith 
1269566063dSJacob Faibussowitsch         PetscCall(MatCreateVecs(pc->pmat, &x, NULL));
1279566063dSJacob Faibussowitsch         PetscCall(VecGetSize(x, &M));
1289566063dSJacob Faibussowitsch         PetscCall(VecGetOwnershipRange(x, &mstart, &mend));
1291b81debcSHong Zhang         mlocal = mend - mstart;
1309566063dSJacob Faibussowitsch         PetscCall(PetscMalloc2(red->psubcomm->n * mlocal, &idx1, red->psubcomm->n * mlocal, &idx2));
1313f457be1SHong Zhang         j = 0;
1321fbd8f88SHong Zhang         for (k = 0; k < red->psubcomm->n; k++) {
1333f457be1SHong Zhang           for (i = mstart; i < mend; i++) {
1343f457be1SHong Zhang             idx1[j]   = i;
135ddc54837SHong Zhang             idx2[j++] = i + M * k;
1363f457be1SHong Zhang           }
1373f457be1SHong Zhang         }
1389566063dSJacob Faibussowitsch         PetscCall(ISCreateGeneral(comm, red->psubcomm->n * mlocal, idx1, PETSC_COPY_VALUES, &is1));
1399566063dSJacob Faibussowitsch         PetscCall(ISCreateGeneral(comm, red->psubcomm->n * mlocal, idx2, PETSC_COPY_VALUES, &is2));
1409566063dSJacob Faibussowitsch         PetscCall(VecScatterCreate(x, is1, red->xdup, is2, &red->scatterin));
1419566063dSJacob Faibussowitsch         PetscCall(ISDestroy(&is1));
1429566063dSJacob Faibussowitsch         PetscCall(ISDestroy(&is2));
1433f457be1SHong Zhang 
1446909748bSHong Zhang         /* Impl below is good for PETSC_SUBCOMM_INTERLACED (no inter-process communication) and PETSC_SUBCOMM_CONTIGUOUS (communication within subcomm) */
1459566063dSJacob Faibussowitsch         PetscCall(ISCreateStride(comm, mlocal, mstart + red->psubcomm->color * M, 1, &is1));
1469566063dSJacob Faibussowitsch         PetscCall(ISCreateStride(comm, mlocal, mstart, 1, &is2));
1479566063dSJacob Faibussowitsch         PetscCall(VecScatterCreate(red->xdup, is1, x, is2, &red->scatterout));
1489566063dSJacob Faibussowitsch         PetscCall(ISDestroy(&is1));
1499566063dSJacob Faibussowitsch         PetscCall(ISDestroy(&is2));
1509566063dSJacob Faibussowitsch         PetscCall(PetscFree2(idx1, idx2));
1519566063dSJacob Faibussowitsch         PetscCall(VecDestroy(&x));
1521b81debcSHong Zhang       }
153ab661555SHong Zhang     } else { /* !red->useparallelmat */
1549566063dSJacob Faibussowitsch       PetscCall(KSPSetOperators(red->ksp, pc->mat, pc->pmat));
1551b81debcSHong Zhang     }
156ab661555SHong Zhang   } else { /* pc->setupcalled */
1574b9ad928SBarry Smith     if (red->useparallelmat) {
158ab661555SHong Zhang       MatReuse reuse;
1591b81debcSHong Zhang       /* grab the parallel matrix and put it into processors of a subcomminicator */
1601b81debcSHong Zhang       /*--------------------------------------------------------------------------*/
161ab661555SHong Zhang       if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1624b9ad928SBarry Smith         /* destroy old matrices */
1639566063dSJacob Faibussowitsch         PetscCall(MatDestroy(&red->pmats));
164ab661555SHong Zhang         reuse = MAT_INITIAL_MATRIX;
1654b9ad928SBarry Smith       } else {
166ab661555SHong Zhang         reuse = MAT_REUSE_MATRIX;
167ab661555SHong Zhang       }
1689566063dSJacob Faibussowitsch       PetscCall(MatCreateRedundantMatrix(pc->pmat, red->psubcomm->n, PetscSubcommChild(red->psubcomm), reuse, &red->pmats));
1699566063dSJacob Faibussowitsch       PetscCall(KSPSetOperators(red->ksp, red->pmats, red->pmats));
170ab661555SHong Zhang     } else { /* !red->useparallelmat */
1719566063dSJacob Faibussowitsch       PetscCall(KSPSetOperators(red->ksp, pc->mat, pc->pmat));
1724b9ad928SBarry Smith     }
173ab661555SHong Zhang   }
1741b81debcSHong Zhang 
1751baa6e33SBarry Smith   if (pc->setfromoptionscalled) PetscCall(KSPSetFromOptions(red->ksp));
1769566063dSJacob Faibussowitsch   PetscCall(KSPSetUp(red->ksp));
1774b9ad928SBarry Smith   PetscFunctionReturn(0);
1784b9ad928SBarry Smith }
1794b9ad928SBarry Smith 
1809371c9d4SSatish Balay static PetscErrorCode PCApply_Redundant(PC pc, Vec x, Vec y) {
1814b9ad928SBarry Smith   PC_Redundant *red = (PC_Redundant *)pc->data;
1823f457be1SHong Zhang   PetscScalar  *array;
1834b9ad928SBarry Smith 
1844b9ad928SBarry Smith   PetscFunctionBegin;
185ddc54837SHong Zhang   if (!red->useparallelmat) {
1869566063dSJacob Faibussowitsch     PetscCall(KSPSolve(red->ksp, x, y));
1879566063dSJacob Faibussowitsch     PetscCall(KSPCheckSolve(red->ksp, pc, y));
188ddc54837SHong Zhang     PetscFunctionReturn(0);
189ddc54837SHong Zhang   }
190ddc54837SHong Zhang 
1913f457be1SHong Zhang   /* scatter x to xdup */
1929566063dSJacob Faibussowitsch   PetscCall(VecScatterBegin(red->scatterin, x, red->xdup, INSERT_VALUES, SCATTER_FORWARD));
1939566063dSJacob Faibussowitsch   PetscCall(VecScatterEnd(red->scatterin, x, red->xdup, INSERT_VALUES, SCATTER_FORWARD));
1943f457be1SHong Zhang 
1953f457be1SHong Zhang   /* place xdup's local array into xsub */
1969566063dSJacob Faibussowitsch   PetscCall(VecGetArray(red->xdup, &array));
1979566063dSJacob Faibussowitsch   PetscCall(VecPlaceArray(red->xsub, (const PetscScalar *)array));
1984b9ad928SBarry Smith 
1994b9ad928SBarry Smith   /* apply preconditioner on each processor */
2009566063dSJacob Faibussowitsch   PetscCall(KSPSolve(red->ksp, red->xsub, red->ysub));
2019566063dSJacob Faibussowitsch   PetscCall(KSPCheckSolve(red->ksp, pc, red->ysub));
2029566063dSJacob Faibussowitsch   PetscCall(VecResetArray(red->xsub));
2039566063dSJacob Faibussowitsch   PetscCall(VecRestoreArray(red->xdup, &array));
2044b9ad928SBarry Smith 
2053f457be1SHong Zhang   /* place ysub's local array into ydup */
2069566063dSJacob Faibussowitsch   PetscCall(VecGetArray(red->ysub, &array));
2079566063dSJacob Faibussowitsch   PetscCall(VecPlaceArray(red->ydup, (const PetscScalar *)array));
2083f457be1SHong Zhang 
2093f457be1SHong Zhang   /* scatter ydup to y */
2109566063dSJacob Faibussowitsch   PetscCall(VecScatterBegin(red->scatterout, red->ydup, y, INSERT_VALUES, SCATTER_FORWARD));
2119566063dSJacob Faibussowitsch   PetscCall(VecScatterEnd(red->scatterout, red->ydup, y, INSERT_VALUES, SCATTER_FORWARD));
2129566063dSJacob Faibussowitsch   PetscCall(VecResetArray(red->ydup));
2139566063dSJacob Faibussowitsch   PetscCall(VecRestoreArray(red->ysub, &array));
2144b9ad928SBarry Smith   PetscFunctionReturn(0);
2154b9ad928SBarry Smith }
2164b9ad928SBarry Smith 
2179371c9d4SSatish Balay static PetscErrorCode PCApplyTranspose_Redundant(PC pc, Vec x, Vec y) {
218d88bfacbSStefano Zampini   PC_Redundant *red = (PC_Redundant *)pc->data;
219d88bfacbSStefano Zampini   PetscScalar  *array;
220d88bfacbSStefano Zampini 
221d88bfacbSStefano Zampini   PetscFunctionBegin;
222d88bfacbSStefano Zampini   if (!red->useparallelmat) {
2239566063dSJacob Faibussowitsch     PetscCall(KSPSolveTranspose(red->ksp, x, y));
2249566063dSJacob Faibussowitsch     PetscCall(KSPCheckSolve(red->ksp, pc, y));
225d88bfacbSStefano Zampini     PetscFunctionReturn(0);
226d88bfacbSStefano Zampini   }
227d88bfacbSStefano Zampini 
228d88bfacbSStefano Zampini   /* scatter x to xdup */
2299566063dSJacob Faibussowitsch   PetscCall(VecScatterBegin(red->scatterin, x, red->xdup, INSERT_VALUES, SCATTER_FORWARD));
2309566063dSJacob Faibussowitsch   PetscCall(VecScatterEnd(red->scatterin, x, red->xdup, INSERT_VALUES, SCATTER_FORWARD));
231d88bfacbSStefano Zampini 
232d88bfacbSStefano Zampini   /* place xdup's local array into xsub */
2339566063dSJacob Faibussowitsch   PetscCall(VecGetArray(red->xdup, &array));
2349566063dSJacob Faibussowitsch   PetscCall(VecPlaceArray(red->xsub, (const PetscScalar *)array));
235d88bfacbSStefano Zampini 
236d88bfacbSStefano Zampini   /* apply preconditioner on each processor */
2379566063dSJacob Faibussowitsch   PetscCall(KSPSolveTranspose(red->ksp, red->xsub, red->ysub));
2389566063dSJacob Faibussowitsch   PetscCall(KSPCheckSolve(red->ksp, pc, red->ysub));
2399566063dSJacob Faibussowitsch   PetscCall(VecResetArray(red->xsub));
2409566063dSJacob Faibussowitsch   PetscCall(VecRestoreArray(red->xdup, &array));
241d88bfacbSStefano Zampini 
242d88bfacbSStefano Zampini   /* place ysub's local array into ydup */
2439566063dSJacob Faibussowitsch   PetscCall(VecGetArray(red->ysub, &array));
2449566063dSJacob Faibussowitsch   PetscCall(VecPlaceArray(red->ydup, (const PetscScalar *)array));
245d88bfacbSStefano Zampini 
246d88bfacbSStefano Zampini   /* scatter ydup to y */
2479566063dSJacob Faibussowitsch   PetscCall(VecScatterBegin(red->scatterout, red->ydup, y, INSERT_VALUES, SCATTER_FORWARD));
2489566063dSJacob Faibussowitsch   PetscCall(VecScatterEnd(red->scatterout, red->ydup, y, INSERT_VALUES, SCATTER_FORWARD));
2499566063dSJacob Faibussowitsch   PetscCall(VecResetArray(red->ydup));
2509566063dSJacob Faibussowitsch   PetscCall(VecRestoreArray(red->ysub, &array));
251d88bfacbSStefano Zampini   PetscFunctionReturn(0);
252d88bfacbSStefano Zampini }
253d88bfacbSStefano Zampini 
2549371c9d4SSatish Balay static PetscErrorCode PCReset_Redundant(PC pc) {
2554b9ad928SBarry Smith   PC_Redundant *red = (PC_Redundant *)pc->data;
2564b9ad928SBarry Smith 
2574b9ad928SBarry Smith   PetscFunctionBegin;
2581b81debcSHong Zhang   if (red->useparallelmat) {
2599566063dSJacob Faibussowitsch     PetscCall(VecScatterDestroy(&red->scatterin));
2609566063dSJacob Faibussowitsch     PetscCall(VecScatterDestroy(&red->scatterout));
2619566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&red->ysub));
2629566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&red->xsub));
2639566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&red->xdup));
2649566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&red->ydup));
2651b81debcSHong Zhang   }
2669566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&red->pmats));
2679566063dSJacob Faibussowitsch   PetscCall(KSPReset(red->ksp));
2681ea5a559SBarry Smith   PetscFunctionReturn(0);
2691ea5a559SBarry Smith }
2701ea5a559SBarry Smith 
2719371c9d4SSatish Balay static PetscErrorCode PCDestroy_Redundant(PC pc) {
2721ea5a559SBarry Smith   PC_Redundant *red = (PC_Redundant *)pc->data;
2731ea5a559SBarry Smith 
2741ea5a559SBarry Smith   PetscFunctionBegin;
2759566063dSJacob Faibussowitsch   PetscCall(PCReset_Redundant(pc));
2769566063dSJacob Faibussowitsch   PetscCall(KSPDestroy(&red->ksp));
2779566063dSJacob Faibussowitsch   PetscCall(PetscSubcommDestroy(&red->psubcomm));
2782e956fe4SStefano Zampini   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantSetScatter_C", NULL));
2792e956fe4SStefano Zampini   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantSetNumber_C", NULL));
2802e956fe4SStefano Zampini   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantGetKSP_C", NULL));
2812e956fe4SStefano Zampini   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantGetOperators_C", NULL));
2822e956fe4SStefano Zampini   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFactorSetShiftType_C", NULL));
2839566063dSJacob Faibussowitsch   PetscCall(PetscFree(pc->data));
2844b9ad928SBarry Smith   PetscFunctionReturn(0);
2854b9ad928SBarry Smith }
2864b9ad928SBarry Smith 
2879371c9d4SSatish Balay static PetscErrorCode PCSetFromOptions_Redundant(PC pc, PetscOptionItems *PetscOptionsObject) {
288a98ce0f4SHong Zhang   PC_Redundant *red = (PC_Redundant *)pc->data;
289a98ce0f4SHong Zhang 
2904b9ad928SBarry Smith   PetscFunctionBegin;
291d0609cedSBarry Smith   PetscOptionsHeadBegin(PetscOptionsObject, "Redundant options");
2929566063dSJacob Faibussowitsch   PetscCall(PetscOptionsInt("-pc_redundant_number", "Number of redundant pc", "PCRedundantSetNumber", red->nsubcomm, &red->nsubcomm, NULL));
293d0609cedSBarry Smith   PetscOptionsHeadEnd();
2944b9ad928SBarry Smith   PetscFunctionReturn(0);
2954b9ad928SBarry Smith }
2964b9ad928SBarry Smith 
2979371c9d4SSatish Balay static PetscErrorCode PCRedundantSetNumber_Redundant(PC pc, PetscInt nreds) {
29809a6bc64SHong Zhang   PC_Redundant *red = (PC_Redundant *)pc->data;
29909a6bc64SHong Zhang 
30009a6bc64SHong Zhang   PetscFunctionBegin;
30109a6bc64SHong Zhang   red->nsubcomm = nreds;
30209a6bc64SHong Zhang   PetscFunctionReturn(0);
30309a6bc64SHong Zhang }
30409a6bc64SHong Zhang 
30509a6bc64SHong Zhang /*@
30609a6bc64SHong Zhang    PCRedundantSetNumber - Sets the number of redundant preconditioner contexts.
30709a6bc64SHong Zhang 
308f1580f4eSBarry Smith    Logically Collective on pc
30909a6bc64SHong Zhang 
31009a6bc64SHong Zhang    Input Parameters:
31109a6bc64SHong Zhang +  pc - the preconditioner context
3129b21d695SBarry Smith -  nredundant - number of redundant preconditioner contexts; for example if you are using 64 MPI processes and
3139b21d695SBarry Smith                               use an nredundant of 4 there will be 4 parallel solves each on 16 = 64/4 processes.
31409a6bc64SHong Zhang 
31509a6bc64SHong Zhang    Level: advanced
31609a6bc64SHong Zhang 
317f1580f4eSBarry Smith .seealso: `PCREDUNDANT`
31809a6bc64SHong Zhang @*/
3199371c9d4SSatish Balay PetscErrorCode PCRedundantSetNumber(PC pc, PetscInt nredundant) {
32009a6bc64SHong Zhang   PetscFunctionBegin;
3210700a824SBarry Smith   PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
32263a3b9bcSJacob Faibussowitsch   PetscCheck(nredundant > 0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "num of redundant pc %" PetscInt_FMT " must be positive", nredundant);
323cac4c232SBarry Smith   PetscTryMethod(pc, "PCRedundantSetNumber_C", (PC, PetscInt), (pc, nredundant));
32409a6bc64SHong Zhang   PetscFunctionReturn(0);
32509a6bc64SHong Zhang }
32609a6bc64SHong Zhang 
3279371c9d4SSatish Balay static PetscErrorCode PCRedundantSetScatter_Redundant(PC pc, VecScatter in, VecScatter out) {
3284b9ad928SBarry Smith   PC_Redundant *red = (PC_Redundant *)pc->data;
3294b9ad928SBarry Smith 
3304b9ad928SBarry Smith   PetscFunctionBegin;
3319566063dSJacob Faibussowitsch   PetscCall(PetscObjectReference((PetscObject)in));
3329566063dSJacob Faibussowitsch   PetscCall(VecScatterDestroy(&red->scatterin));
3332fa5cd67SKarl Rupp 
334c3122656SLisandro Dalcin   red->scatterin = in;
3352fa5cd67SKarl Rupp 
3369566063dSJacob Faibussowitsch   PetscCall(PetscObjectReference((PetscObject)out));
3379566063dSJacob Faibussowitsch   PetscCall(VecScatterDestroy(&red->scatterout));
338c3122656SLisandro Dalcin   red->scatterout = out;
3394b9ad928SBarry Smith   PetscFunctionReturn(0);
3404b9ad928SBarry Smith }
3414b9ad928SBarry Smith 
3424b9ad928SBarry Smith /*@
3434b9ad928SBarry Smith    PCRedundantSetScatter - Sets the scatter used to copy values into the
3444b9ad928SBarry Smith      redundant local solve and the scatter to move them back into the global
3454b9ad928SBarry Smith      vector.
3464b9ad928SBarry Smith 
347f1580f4eSBarry Smith    Logically Collective on pc
3484b9ad928SBarry Smith 
3494b9ad928SBarry Smith    Input Parameters:
3504b9ad928SBarry Smith +  pc - the preconditioner context
3514b9ad928SBarry Smith .  in - the scatter to move the values in
3524b9ad928SBarry Smith -  out - the scatter to move them out
3534b9ad928SBarry Smith 
3544b9ad928SBarry Smith    Level: advanced
3554b9ad928SBarry Smith 
356f1580f4eSBarry Smith .seealso: `PCREDUNDANT`
3574b9ad928SBarry Smith @*/
3589371c9d4SSatish Balay PetscErrorCode PCRedundantSetScatter(PC pc, VecScatter in, VecScatter out) {
3594b9ad928SBarry Smith   PetscFunctionBegin;
3600700a824SBarry Smith   PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
36197929ea7SJunchao Zhang   PetscValidHeaderSpecific(in, PETSCSF_CLASSID, 2);
36297929ea7SJunchao Zhang   PetscValidHeaderSpecific(out, PETSCSF_CLASSID, 3);
363cac4c232SBarry Smith   PetscTryMethod(pc, "PCRedundantSetScatter_C", (PC, VecScatter, VecScatter), (pc, in, out));
3644b9ad928SBarry Smith   PetscFunctionReturn(0);
3654b9ad928SBarry Smith }
3664b9ad928SBarry Smith 
3679371c9d4SSatish Balay static PetscErrorCode PCRedundantGetKSP_Redundant(PC pc, KSP *innerksp) {
3684b9ad928SBarry Smith   PC_Redundant *red = (PC_Redundant *)pc->data;
36975024027SHong Zhang   MPI_Comm      comm, subcomm;
37075024027SHong Zhang   const char   *prefix;
37108cecb0aSPierre Jolivet   PetscBool     issbaij;
3724b9ad928SBarry Smith 
3734b9ad928SBarry Smith   PetscFunctionBegin;
37475024027SHong Zhang   if (!red->psubcomm) {
3759566063dSJacob Faibussowitsch     PetscCall(PCGetOptionsPrefix(pc, &prefix));
376e5acf8a4SHong Zhang 
3779566063dSJacob Faibussowitsch     PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
3789566063dSJacob Faibussowitsch     PetscCall(PetscSubcommCreate(comm, &red->psubcomm));
3799566063dSJacob Faibussowitsch     PetscCall(PetscSubcommSetNumber(red->psubcomm, red->nsubcomm));
3809566063dSJacob Faibussowitsch     PetscCall(PetscSubcommSetType(red->psubcomm, PETSC_SUBCOMM_CONTIGUOUS));
381e5acf8a4SHong Zhang 
3829566063dSJacob Faibussowitsch     PetscCall(PetscSubcommSetOptionsPrefix(red->psubcomm, prefix));
3839566063dSJacob Faibussowitsch     PetscCall(PetscSubcommSetFromOptions(red->psubcomm));
38475024027SHong Zhang 
38575024027SHong Zhang     /* create a new PC that processors in each subcomm have copy of */
38675024027SHong Zhang     subcomm = PetscSubcommChild(red->psubcomm);
38775024027SHong Zhang 
3889566063dSJacob Faibussowitsch     PetscCall(KSPCreate(subcomm, &red->ksp));
3899566063dSJacob Faibussowitsch     PetscCall(KSPSetErrorIfNotConverged(red->ksp, pc->erroriffailure));
3909566063dSJacob Faibussowitsch     PetscCall(PetscObjectIncrementTabLevel((PetscObject)red->ksp, (PetscObject)pc, 1));
3919566063dSJacob Faibussowitsch     PetscCall(KSPSetType(red->ksp, KSPPREONLY));
3929566063dSJacob Faibussowitsch     PetscCall(KSPGetPC(red->ksp, &red->pc));
3939566063dSJacob Faibussowitsch     PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATSEQSBAIJ, &issbaij));
39448a46eb9SPierre Jolivet     if (!issbaij) PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATMPISBAIJ, &issbaij));
39508cecb0aSPierre Jolivet     if (!issbaij) {
3969566063dSJacob Faibussowitsch       PetscCall(PCSetType(red->pc, PCLU));
39708cecb0aSPierre Jolivet     } else {
3989566063dSJacob Faibussowitsch       PetscCall(PCSetType(red->pc, PCCHOLESKY));
39908cecb0aSPierre Jolivet     }
400753b7fb9SBarry Smith     if (red->shifttypeset) {
4019566063dSJacob Faibussowitsch       PetscCall(PCFactorSetShiftType(red->pc, red->shifttype));
402753b7fb9SBarry Smith       red->shifttypeset = PETSC_FALSE;
403753b7fb9SBarry Smith     }
4049566063dSJacob Faibussowitsch     PetscCall(KSPSetOptionsPrefix(red->ksp, prefix));
4059566063dSJacob Faibussowitsch     PetscCall(KSPAppendOptionsPrefix(red->ksp, "redundant_"));
40675024027SHong Zhang   }
40783ab6a24SBarry Smith   *innerksp = red->ksp;
4084b9ad928SBarry Smith   PetscFunctionReturn(0);
4094b9ad928SBarry Smith }
4104b9ad928SBarry Smith 
4114b9ad928SBarry Smith /*@
412f1580f4eSBarry Smith    PCRedundantGetKSP - Gets the less parallel `KSP` created by the redundant `PC`.
4134b9ad928SBarry Smith 
4144b9ad928SBarry Smith    Not Collective
4154b9ad928SBarry Smith 
4164b9ad928SBarry Smith    Input Parameter:
4174b9ad928SBarry Smith .  pc - the preconditioner context
4184b9ad928SBarry Smith 
4194b9ad928SBarry Smith    Output Parameter:
420f1580f4eSBarry Smith .  innerksp - the `KSP` on the smaller set of processes
4214b9ad928SBarry Smith 
4224b9ad928SBarry Smith    Level: advanced
4234b9ad928SBarry Smith 
424f1580f4eSBarry Smith .seealso: `PCREDUNDANT`
4254b9ad928SBarry Smith @*/
4269371c9d4SSatish Balay PetscErrorCode PCRedundantGetKSP(PC pc, KSP *innerksp) {
4274b9ad928SBarry Smith   PetscFunctionBegin;
4280700a824SBarry Smith   PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
42983ab6a24SBarry Smith   PetscValidPointer(innerksp, 2);
430cac4c232SBarry Smith   PetscUseMethod(pc, "PCRedundantGetKSP_C", (PC, KSP *), (pc, innerksp));
4314b9ad928SBarry Smith   PetscFunctionReturn(0);
4324b9ad928SBarry Smith }
4334b9ad928SBarry Smith 
4349371c9d4SSatish Balay static PetscErrorCode PCRedundantGetOperators_Redundant(PC pc, Mat *mat, Mat *pmat) {
4354b9ad928SBarry Smith   PC_Redundant *red = (PC_Redundant *)pc->data;
4364b9ad928SBarry Smith 
4374b9ad928SBarry Smith   PetscFunctionBegin;
438b3804887SHong Zhang   if (mat) *mat = red->pmats;
439b3804887SHong Zhang   if (pmat) *pmat = red->pmats;
4404b9ad928SBarry Smith   PetscFunctionReturn(0);
4414b9ad928SBarry Smith }
4424b9ad928SBarry Smith 
4434b9ad928SBarry Smith /*@
4444b9ad928SBarry Smith    PCRedundantGetOperators - gets the sequential matrix and preconditioner matrix
4454b9ad928SBarry Smith 
4464b9ad928SBarry Smith    Not Collective
4474b9ad928SBarry Smith 
4484b9ad928SBarry Smith    Input Parameter:
4494b9ad928SBarry Smith .  pc - the preconditioner context
4504b9ad928SBarry Smith 
4514b9ad928SBarry Smith    Output Parameters:
4524b9ad928SBarry Smith +  mat - the matrix
4534b9ad928SBarry Smith -  pmat - the (possibly different) preconditioner matrix
4544b9ad928SBarry Smith 
4554b9ad928SBarry Smith    Level: advanced
4564b9ad928SBarry Smith 
457f1580f4eSBarry Smith .seealso: `PCREDUNDANT`
4584b9ad928SBarry Smith @*/
4599371c9d4SSatish Balay PetscErrorCode PCRedundantGetOperators(PC pc, Mat *mat, Mat *pmat) {
4604b9ad928SBarry Smith   PetscFunctionBegin;
4610700a824SBarry Smith   PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
4624482741eSBarry Smith   if (mat) PetscValidPointer(mat, 2);
4634482741eSBarry Smith   if (pmat) PetscValidPointer(pmat, 3);
464cac4c232SBarry Smith   PetscUseMethod(pc, "PCRedundantGetOperators_C", (PC, Mat *, Mat *), (pc, mat, pmat));
4654b9ad928SBarry Smith   PetscFunctionReturn(0);
4664b9ad928SBarry Smith }
4674b9ad928SBarry Smith 
46837a17b4dSBarry Smith /*MC
469f1580f4eSBarry Smith      PCREDUNDANT - Runs a `KSP` solver with preconditioner for the entire problem on subgroups of processors
47037a17b4dSBarry Smith 
471f1580f4eSBarry Smith   Options Database Key:
4729b21d695SBarry Smith .  -pc_redundant_number <n> - number of redundant solves, for example if you are using 64 MPI processes and
4739b21d695SBarry Smith                               use an n of 4 there will be 4 parallel solves each on 16 = 64/4 processes.
47409391456SBarry Smith 
47537a17b4dSBarry Smith    Level: intermediate
47637a17b4dSBarry Smith 
47795452b02SPatrick Sanan    Notes:
478f1580f4eSBarry Smith    Options for the redundant preconditioners can be set using the options database prefix -redundant_
47983ab6a24SBarry Smith 
480f1580f4eSBarry Smith    The default `KSP` is preonly and the default `PC` is `PCLU` or `PCCHOLESKY` if Pmat is of type `MATSBAIJ`.
481753b7fb9SBarry Smith 
482f1580f4eSBarry Smith    `PCFactorSetShiftType()` applied to this `PC` will convey they shift type into the inner `PC` if it is factorization based.
483f1580f4eSBarry Smith 
484f1580f4eSBarry Smith    Developer Note:
485f1580f4eSBarry Smith    `PCSetInitialGuessNonzero()` is not used by this class but likely should be.
4869cfaa89bSBarry Smith 
487db781477SPatrick Sanan .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PCRedundantSetScatter()`,
488f1580f4eSBarry Smith           `PCRedundantGetKSP()`, `PCRedundantGetOperators()`, `PCRedundantSetNumber()`, `PCREDISTRIBUTE`
48937a17b4dSBarry Smith M*/
49037a17b4dSBarry Smith 
4919371c9d4SSatish Balay PETSC_EXTERN PetscErrorCode PCCreate_Redundant(PC pc) {
4924b9ad928SBarry Smith   PC_Redundant *red;
49369db28dcSHong Zhang   PetscMPIInt   size;
4943f457be1SHong Zhang 
4954b9ad928SBarry Smith   PetscFunctionBegin;
496*4dfa11a4SJacob Faibussowitsch   PetscCall(PetscNew(&red));
4979566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
4982fa5cd67SKarl Rupp 
49969db28dcSHong Zhang   red->nsubcomm       = size;
5004b9ad928SBarry Smith   red->useparallelmat = PETSC_TRUE;
5011fbd8f88SHong Zhang   pc->data            = (void *)red;
5024b9ad928SBarry Smith 
5034b9ad928SBarry Smith   pc->ops->apply          = PCApply_Redundant;
504d88bfacbSStefano Zampini   pc->ops->applytranspose = PCApplyTranspose_Redundant;
5054b9ad928SBarry Smith   pc->ops->setup          = PCSetUp_Redundant;
5064b9ad928SBarry Smith   pc->ops->destroy        = PCDestroy_Redundant;
5071ea5a559SBarry Smith   pc->ops->reset          = PCReset_Redundant;
5084b9ad928SBarry Smith   pc->ops->setfromoptions = PCSetFromOptions_Redundant;
5094b9ad928SBarry Smith   pc->ops->view           = PCView_Redundant;
5102fa5cd67SKarl Rupp 
5119566063dSJacob Faibussowitsch   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantSetScatter_C", PCRedundantSetScatter_Redundant));
5129566063dSJacob Faibussowitsch   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantSetNumber_C", PCRedundantSetNumber_Redundant));
5139566063dSJacob Faibussowitsch   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantGetKSP_C", PCRedundantGetKSP_Redundant));
5149566063dSJacob Faibussowitsch   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantGetOperators_C", PCRedundantGetOperators_Redundant));
5159566063dSJacob Faibussowitsch   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFactorSetShiftType_C", PCFactorSetShiftType_Redundant));
5164b9ad928SBarry Smith   PetscFunctionReturn(0);
5174b9ad928SBarry Smith }
518