xref: /petsc/src/ksp/pc/impls/redundant/redundant.c (revision f1580f4e3ce5d5b2393648fd039d0d41b440385d)
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 
308*f1580f4eSBarry 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 
317*f1580f4eSBarry 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 
347*f1580f4eSBarry 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 
356*f1580f4eSBarry 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));
3849566063dSJacob Faibussowitsch     PetscCall(PetscLogObjectMemory((PetscObject)pc, sizeof(PetscSubcomm)));
38575024027SHong Zhang 
38675024027SHong Zhang     /* create a new PC that processors in each subcomm have copy of */
38775024027SHong Zhang     subcomm = PetscSubcommChild(red->psubcomm);
38875024027SHong Zhang 
3899566063dSJacob Faibussowitsch     PetscCall(KSPCreate(subcomm, &red->ksp));
3909566063dSJacob Faibussowitsch     PetscCall(KSPSetErrorIfNotConverged(red->ksp, pc->erroriffailure));
3919566063dSJacob Faibussowitsch     PetscCall(PetscObjectIncrementTabLevel((PetscObject)red->ksp, (PetscObject)pc, 1));
3929566063dSJacob Faibussowitsch     PetscCall(PetscLogObjectParent((PetscObject)pc, (PetscObject)red->ksp));
3939566063dSJacob Faibussowitsch     PetscCall(KSPSetType(red->ksp, KSPPREONLY));
3949566063dSJacob Faibussowitsch     PetscCall(KSPGetPC(red->ksp, &red->pc));
3959566063dSJacob Faibussowitsch     PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATSEQSBAIJ, &issbaij));
39648a46eb9SPierre Jolivet     if (!issbaij) PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATMPISBAIJ, &issbaij));
39708cecb0aSPierre Jolivet     if (!issbaij) {
3989566063dSJacob Faibussowitsch       PetscCall(PCSetType(red->pc, PCLU));
39908cecb0aSPierre Jolivet     } else {
4009566063dSJacob Faibussowitsch       PetscCall(PCSetType(red->pc, PCCHOLESKY));
40108cecb0aSPierre Jolivet     }
402753b7fb9SBarry Smith     if (red->shifttypeset) {
4039566063dSJacob Faibussowitsch       PetscCall(PCFactorSetShiftType(red->pc, red->shifttype));
404753b7fb9SBarry Smith       red->shifttypeset = PETSC_FALSE;
405753b7fb9SBarry Smith     }
4069566063dSJacob Faibussowitsch     PetscCall(KSPSetOptionsPrefix(red->ksp, prefix));
4079566063dSJacob Faibussowitsch     PetscCall(KSPAppendOptionsPrefix(red->ksp, "redundant_"));
40875024027SHong Zhang   }
40983ab6a24SBarry Smith   *innerksp = red->ksp;
4104b9ad928SBarry Smith   PetscFunctionReturn(0);
4114b9ad928SBarry Smith }
4124b9ad928SBarry Smith 
4134b9ad928SBarry Smith /*@
414*f1580f4eSBarry Smith    PCRedundantGetKSP - Gets the less parallel `KSP` created by the redundant `PC`.
4154b9ad928SBarry Smith 
4164b9ad928SBarry Smith    Not Collective
4174b9ad928SBarry Smith 
4184b9ad928SBarry Smith    Input Parameter:
4194b9ad928SBarry Smith .  pc - the preconditioner context
4204b9ad928SBarry Smith 
4214b9ad928SBarry Smith    Output Parameter:
422*f1580f4eSBarry Smith .  innerksp - the `KSP` on the smaller set of processes
4234b9ad928SBarry Smith 
4244b9ad928SBarry Smith    Level: advanced
4254b9ad928SBarry Smith 
426*f1580f4eSBarry Smith .seealso: `PCREDUNDANT`
4274b9ad928SBarry Smith @*/
4289371c9d4SSatish Balay PetscErrorCode PCRedundantGetKSP(PC pc, KSP *innerksp) {
4294b9ad928SBarry Smith   PetscFunctionBegin;
4300700a824SBarry Smith   PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
43183ab6a24SBarry Smith   PetscValidPointer(innerksp, 2);
432cac4c232SBarry Smith   PetscUseMethod(pc, "PCRedundantGetKSP_C", (PC, KSP *), (pc, innerksp));
4334b9ad928SBarry Smith   PetscFunctionReturn(0);
4344b9ad928SBarry Smith }
4354b9ad928SBarry Smith 
4369371c9d4SSatish Balay static PetscErrorCode PCRedundantGetOperators_Redundant(PC pc, Mat *mat, Mat *pmat) {
4374b9ad928SBarry Smith   PC_Redundant *red = (PC_Redundant *)pc->data;
4384b9ad928SBarry Smith 
4394b9ad928SBarry Smith   PetscFunctionBegin;
440b3804887SHong Zhang   if (mat) *mat = red->pmats;
441b3804887SHong Zhang   if (pmat) *pmat = red->pmats;
4424b9ad928SBarry Smith   PetscFunctionReturn(0);
4434b9ad928SBarry Smith }
4444b9ad928SBarry Smith 
4454b9ad928SBarry Smith /*@
4464b9ad928SBarry Smith    PCRedundantGetOperators - gets the sequential matrix and preconditioner matrix
4474b9ad928SBarry Smith 
4484b9ad928SBarry Smith    Not Collective
4494b9ad928SBarry Smith 
4504b9ad928SBarry Smith    Input Parameter:
4514b9ad928SBarry Smith .  pc - the preconditioner context
4524b9ad928SBarry Smith 
4534b9ad928SBarry Smith    Output Parameters:
4544b9ad928SBarry Smith +  mat - the matrix
4554b9ad928SBarry Smith -  pmat - the (possibly different) preconditioner matrix
4564b9ad928SBarry Smith 
4574b9ad928SBarry Smith    Level: advanced
4584b9ad928SBarry Smith 
459*f1580f4eSBarry Smith .seealso: `PCREDUNDANT`
4604b9ad928SBarry Smith @*/
4619371c9d4SSatish Balay PetscErrorCode PCRedundantGetOperators(PC pc, Mat *mat, Mat *pmat) {
4624b9ad928SBarry Smith   PetscFunctionBegin;
4630700a824SBarry Smith   PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
4644482741eSBarry Smith   if (mat) PetscValidPointer(mat, 2);
4654482741eSBarry Smith   if (pmat) PetscValidPointer(pmat, 3);
466cac4c232SBarry Smith   PetscUseMethod(pc, "PCRedundantGetOperators_C", (PC, Mat *, Mat *), (pc, mat, pmat));
4674b9ad928SBarry Smith   PetscFunctionReturn(0);
4684b9ad928SBarry Smith }
4694b9ad928SBarry Smith 
47037a17b4dSBarry Smith /*MC
471*f1580f4eSBarry Smith      PCREDUNDANT - Runs a `KSP` solver with preconditioner for the entire problem on subgroups of processors
47237a17b4dSBarry Smith 
473*f1580f4eSBarry Smith   Options Database Key:
4749b21d695SBarry Smith .  -pc_redundant_number <n> - number of redundant solves, for example if you are using 64 MPI processes and
4759b21d695SBarry Smith                               use an n of 4 there will be 4 parallel solves each on 16 = 64/4 processes.
47609391456SBarry Smith 
47737a17b4dSBarry Smith    Level: intermediate
47837a17b4dSBarry Smith 
47995452b02SPatrick Sanan    Notes:
480*f1580f4eSBarry Smith    Options for the redundant preconditioners can be set using the options database prefix -redundant_
48183ab6a24SBarry Smith 
482*f1580f4eSBarry Smith    The default `KSP` is preonly and the default `PC` is `PCLU` or `PCCHOLESKY` if Pmat is of type `MATSBAIJ`.
483753b7fb9SBarry Smith 
484*f1580f4eSBarry Smith    `PCFactorSetShiftType()` applied to this `PC` will convey they shift type into the inner `PC` if it is factorization based.
485*f1580f4eSBarry Smith 
486*f1580f4eSBarry Smith    Developer Note:
487*f1580f4eSBarry Smith    `PCSetInitialGuessNonzero()` is not used by this class but likely should be.
4889cfaa89bSBarry Smith 
489db781477SPatrick Sanan .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PCRedundantSetScatter()`,
490*f1580f4eSBarry Smith           `PCRedundantGetKSP()`, `PCRedundantGetOperators()`, `PCRedundantSetNumber()`, `PCREDISTRIBUTE`
49137a17b4dSBarry Smith M*/
49237a17b4dSBarry Smith 
4939371c9d4SSatish Balay PETSC_EXTERN PetscErrorCode PCCreate_Redundant(PC pc) {
4944b9ad928SBarry Smith   PC_Redundant *red;
49569db28dcSHong Zhang   PetscMPIInt   size;
4963f457be1SHong Zhang 
4974b9ad928SBarry Smith   PetscFunctionBegin;
4989566063dSJacob Faibussowitsch   PetscCall(PetscNewLog(pc, &red));
4999566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
5002fa5cd67SKarl Rupp 
50169db28dcSHong Zhang   red->nsubcomm       = size;
5024b9ad928SBarry Smith   red->useparallelmat = PETSC_TRUE;
5031fbd8f88SHong Zhang   pc->data            = (void *)red;
5044b9ad928SBarry Smith 
5054b9ad928SBarry Smith   pc->ops->apply          = PCApply_Redundant;
506d88bfacbSStefano Zampini   pc->ops->applytranspose = PCApplyTranspose_Redundant;
5074b9ad928SBarry Smith   pc->ops->setup          = PCSetUp_Redundant;
5084b9ad928SBarry Smith   pc->ops->destroy        = PCDestroy_Redundant;
5091ea5a559SBarry Smith   pc->ops->reset          = PCReset_Redundant;
5104b9ad928SBarry Smith   pc->ops->setfromoptions = PCSetFromOptions_Redundant;
5114b9ad928SBarry Smith   pc->ops->view           = PCView_Redundant;
5122fa5cd67SKarl Rupp 
5139566063dSJacob Faibussowitsch   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantSetScatter_C", PCRedundantSetScatter_Redundant));
5149566063dSJacob Faibussowitsch   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantSetNumber_C", PCRedundantSetNumber_Redundant));
5159566063dSJacob Faibussowitsch   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantGetKSP_C", PCRedundantGetKSP_Redundant));
5169566063dSJacob Faibussowitsch   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantGetOperators_C", PCRedundantGetOperators_Redundant));
5179566063dSJacob Faibussowitsch   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFactorSetShiftType_C", PCFactorSetShiftType_Redundant));
5184b9ad928SBarry Smith   PetscFunctionReturn(0);
5194b9ad928SBarry Smith }
520