xref: /petsc/src/ksp/pc/impls/redundant/redundant.c (revision 9767e1d8b281455ed5b7d022ec0a9ffdd7f4cae7)
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 
22d71ae5a4SJacob Faibussowitsch PetscErrorCode PCFactorSetShiftType_Redundant(PC pc, MatFactorShiftType shifttype)
23d71ae5a4SJacob Faibussowitsch {
24753b7fb9SBarry Smith   PC_Redundant *red = (PC_Redundant *)pc->data;
25753b7fb9SBarry Smith 
26753b7fb9SBarry Smith   PetscFunctionBegin;
27753b7fb9SBarry Smith   if (red->ksp) {
28753b7fb9SBarry Smith     PC pc;
299566063dSJacob Faibussowitsch     PetscCall(KSPGetPC(red->ksp, &pc));
309566063dSJacob Faibussowitsch     PetscCall(PCFactorSetShiftType(pc, shifttype));
31753b7fb9SBarry Smith   } else {
32753b7fb9SBarry Smith     red->shifttypeset = PETSC_TRUE;
33753b7fb9SBarry Smith     red->shifttype    = shifttype;
34753b7fb9SBarry Smith   }
353ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
36753b7fb9SBarry Smith }
37753b7fb9SBarry Smith 
38d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCView_Redundant(PC pc, PetscViewer viewer)
39d71ae5a4SJacob Faibussowitsch {
404b9ad928SBarry Smith   PC_Redundant *red = (PC_Redundant *)pc->data;
41ace3abfcSBarry Smith   PetscBool     iascii, isstring;
4203ccd0b4SBarry Smith   PetscViewer   subviewer;
434b9ad928SBarry Smith 
444b9ad928SBarry Smith   PetscFunctionBegin;
459566063dSJacob Faibussowitsch   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
469566063dSJacob Faibussowitsch   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
4732077d6dSBarry Smith   if (iascii) {
4803ccd0b4SBarry Smith     if (!red->psubcomm) {
499566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(viewer, "  Not yet setup\n"));
5003ccd0b4SBarry Smith     } else {
5163a3b9bcSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(viewer, "  First (color=0) of %" PetscInt_FMT " PCs follows\n", red->nsubcomm));
529566063dSJacob Faibussowitsch       PetscCall(PetscViewerGetSubViewer(viewer, ((PetscObject)red->pc)->comm, &subviewer));
53f5dd71faSBarry Smith       if (!red->psubcomm->color) { /* only view first redundant pc */
549566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPushTab(subviewer));
559566063dSJacob Faibussowitsch         PetscCall(KSPView(red->ksp, subviewer));
569566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPopTab(subviewer));
574b9ad928SBarry Smith       }
589566063dSJacob Faibussowitsch       PetscCall(PetscViewerRestoreSubViewer(viewer, ((PetscObject)red->pc)->comm, &subviewer));
594b9ad928SBarry Smith     }
6003ccd0b4SBarry Smith   } else if (isstring) {
619566063dSJacob Faibussowitsch     PetscCall(PetscViewerStringSPrintf(viewer, " Redundant solver preconditioner"));
624b9ad928SBarry Smith   }
633ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
644b9ad928SBarry Smith }
654b9ad928SBarry Smith 
6619b3b6edSHong Zhang #include <../src/mat/impls/aij/mpi/mpiaij.h>
67d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCSetUp_Redundant(PC pc)
68d71ae5a4SJacob Faibussowitsch {
694b9ad928SBarry Smith   PC_Redundant *red = (PC_Redundant *)pc->data;
701b81debcSHong Zhang   PetscInt      mstart, mend, mlocal, M;
7113f74950SBarry Smith   PetscMPIInt   size;
72ce94432eSBarry Smith   MPI_Comm      comm, subcomm;
73ddc54837SHong Zhang   Vec           x;
743f457be1SHong Zhang 
754b9ad928SBarry Smith   PetscFunctionBegin;
769566063dSJacob Faibussowitsch   PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
77ddc54837SHong Zhang 
78ddc54837SHong Zhang   /* if pmatrix set by user is sequential then we do not need to gather the parallel matrix */
799566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_size(comm, &size));
80ddc54837SHong Zhang   if (size == 1) red->useparallelmat = PETSC_FALSE;
811fbd8f88SHong Zhang 
824b9ad928SBarry Smith   if (!pc->setupcalled) {
831b81debcSHong Zhang     PetscInt mloc_sub;
8475024027SHong Zhang     if (!red->psubcomm) { /* create red->psubcomm, new ksp and pc over subcomm */
8575024027SHong Zhang       KSP ksp;
869566063dSJacob Faibussowitsch       PetscCall(PCRedundantGetKSP(pc, &ksp));
871b81debcSHong Zhang     }
8875024027SHong Zhang     subcomm = PetscSubcommChild(red->psubcomm);
891fbd8f88SHong Zhang 
901b81debcSHong Zhang     if (red->useparallelmat) {
911b81debcSHong Zhang       /* grab the parallel matrix and put it into processors of a subcomminicator */
929566063dSJacob Faibussowitsch       PetscCall(MatCreateRedundantMatrix(pc->pmat, red->psubcomm->n, subcomm, MAT_INITIAL_MATRIX, &red->pmats));
93b85f2e9bSHong Zhang 
949566063dSJacob Faibussowitsch       PetscCallMPI(MPI_Comm_size(subcomm, &size));
95b85f2e9bSHong Zhang       if (size > 1) {
9608cecb0aSPierre Jolivet         PetscBool foundpack, issbaij;
979566063dSJacob Faibussowitsch         PetscCall(PetscObjectTypeCompare((PetscObject)red->pmats, MATMPISBAIJ, &issbaij));
9808cecb0aSPierre Jolivet         if (!issbaij) {
999566063dSJacob Faibussowitsch           PetscCall(MatGetFactorAvailable(red->pmats, NULL, MAT_FACTOR_LU, &foundpack));
10008cecb0aSPierre Jolivet         } else {
1019566063dSJacob Faibussowitsch           PetscCall(MatGetFactorAvailable(red->pmats, NULL, MAT_FACTOR_CHOLESKY, &foundpack));
10208cecb0aSPierre Jolivet         }
103b85f2e9bSHong Zhang         if (!foundpack) { /* reset default ksp and pc */
1049566063dSJacob Faibussowitsch           PetscCall(KSPSetType(red->ksp, KSPGMRES));
1059566063dSJacob Faibussowitsch           PetscCall(PCSetType(red->pc, PCBJACOBI));
106c1619fb6SBarry Smith         } else {
1079566063dSJacob Faibussowitsch           PetscCall(PCFactorSetMatSolverType(red->pc, NULL));
108b85f2e9bSHong Zhang         }
109b85f2e9bSHong Zhang       }
110b85f2e9bSHong Zhang 
1119566063dSJacob Faibussowitsch       PetscCall(KSPSetOperators(red->ksp, red->pmats, red->pmats));
1124b9ad928SBarry Smith 
1131b81debcSHong Zhang       /* get working vectors xsub and ysub */
1149566063dSJacob Faibussowitsch       PetscCall(MatCreateVecs(red->pmats, &red->xsub, &red->ysub));
1152fa5cd67SKarl Rupp 
1168b96b0d2SHong Zhang       /* create working vectors xdup and ydup.
1178b96b0d2SHong Zhang        xdup concatenates all xsub's contigously to form a mpi vector over dupcomm  (see PetscSubcommCreate_interlaced())
1188b96b0d2SHong Zhang        ydup concatenates all ysub and has empty local arrays because ysub's arrays will be place into it.
119ce94432eSBarry Smith        Note: we use communicator dupcomm, not PetscObjectComm((PetscObject)pc)! */
1209566063dSJacob Faibussowitsch       PetscCall(MatGetLocalSize(red->pmats, &mloc_sub, NULL));
1219566063dSJacob Faibussowitsch       PetscCall(VecCreateMPI(PetscSubcommContiguousParent(red->psubcomm), mloc_sub, PETSC_DECIDE, &red->xdup));
1229566063dSJacob Faibussowitsch       PetscCall(VecCreateMPIWithArray(PetscSubcommContiguousParent(red->psubcomm), 1, mloc_sub, PETSC_DECIDE, NULL, &red->ydup));
1233f457be1SHong Zhang 
124f68be91cSHong Zhang       /* create vecscatters */
125f68be91cSHong Zhang       if (!red->scatterin) { /* efficiency of scatterin is independent from psubcomm_type! */
1263f457be1SHong Zhang         IS        is1, is2;
1273f457be1SHong Zhang         PetscInt *idx1, *idx2, i, j, k;
12845fc02eaSBarry Smith 
1299566063dSJacob Faibussowitsch         PetscCall(MatCreateVecs(pc->pmat, &x, NULL));
1309566063dSJacob Faibussowitsch         PetscCall(VecGetSize(x, &M));
1319566063dSJacob Faibussowitsch         PetscCall(VecGetOwnershipRange(x, &mstart, &mend));
1321b81debcSHong Zhang         mlocal = mend - mstart;
1339566063dSJacob Faibussowitsch         PetscCall(PetscMalloc2(red->psubcomm->n * mlocal, &idx1, red->psubcomm->n * mlocal, &idx2));
1343f457be1SHong Zhang         j = 0;
1351fbd8f88SHong Zhang         for (k = 0; k < red->psubcomm->n; k++) {
1363f457be1SHong Zhang           for (i = mstart; i < mend; i++) {
1373f457be1SHong Zhang             idx1[j]   = i;
138ddc54837SHong Zhang             idx2[j++] = i + M * k;
1393f457be1SHong Zhang           }
1403f457be1SHong Zhang         }
1419566063dSJacob Faibussowitsch         PetscCall(ISCreateGeneral(comm, red->psubcomm->n * mlocal, idx1, PETSC_COPY_VALUES, &is1));
1429566063dSJacob Faibussowitsch         PetscCall(ISCreateGeneral(comm, red->psubcomm->n * mlocal, idx2, PETSC_COPY_VALUES, &is2));
1439566063dSJacob Faibussowitsch         PetscCall(VecScatterCreate(x, is1, red->xdup, is2, &red->scatterin));
1449566063dSJacob Faibussowitsch         PetscCall(ISDestroy(&is1));
1459566063dSJacob Faibussowitsch         PetscCall(ISDestroy(&is2));
1463f457be1SHong Zhang 
1476909748bSHong Zhang         /* Impl below is good for PETSC_SUBCOMM_INTERLACED (no inter-process communication) and PETSC_SUBCOMM_CONTIGUOUS (communication within subcomm) */
1489566063dSJacob Faibussowitsch         PetscCall(ISCreateStride(comm, mlocal, mstart + red->psubcomm->color * M, 1, &is1));
1499566063dSJacob Faibussowitsch         PetscCall(ISCreateStride(comm, mlocal, mstart, 1, &is2));
1509566063dSJacob Faibussowitsch         PetscCall(VecScatterCreate(red->xdup, is1, x, is2, &red->scatterout));
1519566063dSJacob Faibussowitsch         PetscCall(ISDestroy(&is1));
1529566063dSJacob Faibussowitsch         PetscCall(ISDestroy(&is2));
1539566063dSJacob Faibussowitsch         PetscCall(PetscFree2(idx1, idx2));
1549566063dSJacob Faibussowitsch         PetscCall(VecDestroy(&x));
1551b81debcSHong Zhang       }
156ab661555SHong Zhang     } else { /* !red->useparallelmat */
1579566063dSJacob Faibussowitsch       PetscCall(KSPSetOperators(red->ksp, pc->mat, pc->pmat));
1581b81debcSHong Zhang     }
159ab661555SHong Zhang   } else { /* pc->setupcalled */
1604b9ad928SBarry Smith     if (red->useparallelmat) {
161ab661555SHong Zhang       MatReuse reuse;
1621b81debcSHong Zhang       /* grab the parallel matrix and put it into processors of a subcomminicator */
1631b81debcSHong Zhang       /*--------------------------------------------------------------------------*/
164ab661555SHong Zhang       if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1654b9ad928SBarry Smith         /* destroy old matrices */
1669566063dSJacob Faibussowitsch         PetscCall(MatDestroy(&red->pmats));
167ab661555SHong Zhang         reuse = MAT_INITIAL_MATRIX;
1684b9ad928SBarry Smith       } else {
169ab661555SHong Zhang         reuse = MAT_REUSE_MATRIX;
170ab661555SHong Zhang       }
1719566063dSJacob Faibussowitsch       PetscCall(MatCreateRedundantMatrix(pc->pmat, red->psubcomm->n, PetscSubcommChild(red->psubcomm), reuse, &red->pmats));
1729566063dSJacob Faibussowitsch       PetscCall(KSPSetOperators(red->ksp, red->pmats, red->pmats));
173ab661555SHong Zhang     } else { /* !red->useparallelmat */
1749566063dSJacob Faibussowitsch       PetscCall(KSPSetOperators(red->ksp, pc->mat, pc->pmat));
1754b9ad928SBarry Smith     }
176ab661555SHong Zhang   }
1771b81debcSHong Zhang 
1781baa6e33SBarry Smith   if (pc->setfromoptionscalled) PetscCall(KSPSetFromOptions(red->ksp));
1799566063dSJacob Faibussowitsch   PetscCall(KSPSetUp(red->ksp));
180*9767e1d8SStefano Zampini 
181*9767e1d8SStefano Zampini   /* Detect failure */
182*9767e1d8SStefano Zampini   KSPConvergedReason redreason;
183*9767e1d8SStefano Zampini   PetscCall(KSPGetConvergedReason(red->ksp, &redreason));
184*9767e1d8SStefano Zampini   if (redreason) pc->failedreason = PC_SUBPC_ERROR;
1853ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
1864b9ad928SBarry Smith }
1874b9ad928SBarry Smith 
188d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCApply_Redundant(PC pc, Vec x, Vec y)
189d71ae5a4SJacob Faibussowitsch {
1904b9ad928SBarry Smith   PC_Redundant *red = (PC_Redundant *)pc->data;
1913f457be1SHong Zhang   PetscScalar  *array;
1924b9ad928SBarry Smith 
1934b9ad928SBarry Smith   PetscFunctionBegin;
194ddc54837SHong Zhang   if (!red->useparallelmat) {
1959566063dSJacob Faibussowitsch     PetscCall(KSPSolve(red->ksp, x, y));
1969566063dSJacob Faibussowitsch     PetscCall(KSPCheckSolve(red->ksp, pc, y));
1973ba16761SJacob Faibussowitsch     PetscFunctionReturn(PETSC_SUCCESS);
198ddc54837SHong Zhang   }
199ddc54837SHong Zhang 
2003f457be1SHong Zhang   /* scatter x to xdup */
2019566063dSJacob Faibussowitsch   PetscCall(VecScatterBegin(red->scatterin, x, red->xdup, INSERT_VALUES, SCATTER_FORWARD));
2029566063dSJacob Faibussowitsch   PetscCall(VecScatterEnd(red->scatterin, x, red->xdup, INSERT_VALUES, SCATTER_FORWARD));
2033f457be1SHong Zhang 
2043f457be1SHong Zhang   /* place xdup's local array into xsub */
2059566063dSJacob Faibussowitsch   PetscCall(VecGetArray(red->xdup, &array));
2069566063dSJacob Faibussowitsch   PetscCall(VecPlaceArray(red->xsub, (const PetscScalar *)array));
2074b9ad928SBarry Smith 
2084b9ad928SBarry Smith   /* apply preconditioner on each processor */
2099566063dSJacob Faibussowitsch   PetscCall(KSPSolve(red->ksp, red->xsub, red->ysub));
2109566063dSJacob Faibussowitsch   PetscCall(KSPCheckSolve(red->ksp, pc, red->ysub));
2119566063dSJacob Faibussowitsch   PetscCall(VecResetArray(red->xsub));
2129566063dSJacob Faibussowitsch   PetscCall(VecRestoreArray(red->xdup, &array));
2134b9ad928SBarry Smith 
2143f457be1SHong Zhang   /* place ysub's local array into ydup */
2159566063dSJacob Faibussowitsch   PetscCall(VecGetArray(red->ysub, &array));
2169566063dSJacob Faibussowitsch   PetscCall(VecPlaceArray(red->ydup, (const PetscScalar *)array));
2173f457be1SHong Zhang 
2183f457be1SHong Zhang   /* scatter ydup to y */
2199566063dSJacob Faibussowitsch   PetscCall(VecScatterBegin(red->scatterout, red->ydup, y, INSERT_VALUES, SCATTER_FORWARD));
2209566063dSJacob Faibussowitsch   PetscCall(VecScatterEnd(red->scatterout, red->ydup, y, INSERT_VALUES, SCATTER_FORWARD));
2219566063dSJacob Faibussowitsch   PetscCall(VecResetArray(red->ydup));
2229566063dSJacob Faibussowitsch   PetscCall(VecRestoreArray(red->ysub, &array));
2233ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
2244b9ad928SBarry Smith }
2254b9ad928SBarry Smith 
226d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCApplyTranspose_Redundant(PC pc, Vec x, Vec y)
227d71ae5a4SJacob Faibussowitsch {
228d88bfacbSStefano Zampini   PC_Redundant *red = (PC_Redundant *)pc->data;
229d88bfacbSStefano Zampini   PetscScalar  *array;
230d88bfacbSStefano Zampini 
231d88bfacbSStefano Zampini   PetscFunctionBegin;
232d88bfacbSStefano Zampini   if (!red->useparallelmat) {
2339566063dSJacob Faibussowitsch     PetscCall(KSPSolveTranspose(red->ksp, x, y));
2349566063dSJacob Faibussowitsch     PetscCall(KSPCheckSolve(red->ksp, pc, y));
2353ba16761SJacob Faibussowitsch     PetscFunctionReturn(PETSC_SUCCESS);
236d88bfacbSStefano Zampini   }
237d88bfacbSStefano Zampini 
238d88bfacbSStefano Zampini   /* scatter x to xdup */
2399566063dSJacob Faibussowitsch   PetscCall(VecScatterBegin(red->scatterin, x, red->xdup, INSERT_VALUES, SCATTER_FORWARD));
2409566063dSJacob Faibussowitsch   PetscCall(VecScatterEnd(red->scatterin, x, red->xdup, INSERT_VALUES, SCATTER_FORWARD));
241d88bfacbSStefano Zampini 
242d88bfacbSStefano Zampini   /* place xdup's local array into xsub */
2439566063dSJacob Faibussowitsch   PetscCall(VecGetArray(red->xdup, &array));
2449566063dSJacob Faibussowitsch   PetscCall(VecPlaceArray(red->xsub, (const PetscScalar *)array));
245d88bfacbSStefano Zampini 
246d88bfacbSStefano Zampini   /* apply preconditioner on each processor */
2479566063dSJacob Faibussowitsch   PetscCall(KSPSolveTranspose(red->ksp, red->xsub, red->ysub));
2489566063dSJacob Faibussowitsch   PetscCall(KSPCheckSolve(red->ksp, pc, red->ysub));
2499566063dSJacob Faibussowitsch   PetscCall(VecResetArray(red->xsub));
2509566063dSJacob Faibussowitsch   PetscCall(VecRestoreArray(red->xdup, &array));
251d88bfacbSStefano Zampini 
252d88bfacbSStefano Zampini   /* place ysub's local array into ydup */
2539566063dSJacob Faibussowitsch   PetscCall(VecGetArray(red->ysub, &array));
2549566063dSJacob Faibussowitsch   PetscCall(VecPlaceArray(red->ydup, (const PetscScalar *)array));
255d88bfacbSStefano Zampini 
256d88bfacbSStefano Zampini   /* scatter ydup to y */
2579566063dSJacob Faibussowitsch   PetscCall(VecScatterBegin(red->scatterout, red->ydup, y, INSERT_VALUES, SCATTER_FORWARD));
2589566063dSJacob Faibussowitsch   PetscCall(VecScatterEnd(red->scatterout, red->ydup, y, INSERT_VALUES, SCATTER_FORWARD));
2599566063dSJacob Faibussowitsch   PetscCall(VecResetArray(red->ydup));
2609566063dSJacob Faibussowitsch   PetscCall(VecRestoreArray(red->ysub, &array));
2613ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
262d88bfacbSStefano Zampini }
263d88bfacbSStefano Zampini 
264d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCReset_Redundant(PC pc)
265d71ae5a4SJacob Faibussowitsch {
2664b9ad928SBarry Smith   PC_Redundant *red = (PC_Redundant *)pc->data;
2674b9ad928SBarry Smith 
2684b9ad928SBarry Smith   PetscFunctionBegin;
2691b81debcSHong Zhang   if (red->useparallelmat) {
2709566063dSJacob Faibussowitsch     PetscCall(VecScatterDestroy(&red->scatterin));
2719566063dSJacob Faibussowitsch     PetscCall(VecScatterDestroy(&red->scatterout));
2729566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&red->ysub));
2739566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&red->xsub));
2749566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&red->xdup));
2759566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&red->ydup));
2761b81debcSHong Zhang   }
2779566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&red->pmats));
2789566063dSJacob Faibussowitsch   PetscCall(KSPReset(red->ksp));
2793ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
2801ea5a559SBarry Smith }
2811ea5a559SBarry Smith 
282d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCDestroy_Redundant(PC pc)
283d71ae5a4SJacob Faibussowitsch {
2841ea5a559SBarry Smith   PC_Redundant *red = (PC_Redundant *)pc->data;
2851ea5a559SBarry Smith 
2861ea5a559SBarry Smith   PetscFunctionBegin;
2879566063dSJacob Faibussowitsch   PetscCall(PCReset_Redundant(pc));
2889566063dSJacob Faibussowitsch   PetscCall(KSPDestroy(&red->ksp));
2899566063dSJacob Faibussowitsch   PetscCall(PetscSubcommDestroy(&red->psubcomm));
2902e956fe4SStefano Zampini   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantSetScatter_C", NULL));
2912e956fe4SStefano Zampini   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantSetNumber_C", NULL));
2922e956fe4SStefano Zampini   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantGetKSP_C", NULL));
2932e956fe4SStefano Zampini   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantGetOperators_C", NULL));
2942e956fe4SStefano Zampini   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFactorSetShiftType_C", NULL));
2959566063dSJacob Faibussowitsch   PetscCall(PetscFree(pc->data));
2963ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
2974b9ad928SBarry Smith }
2984b9ad928SBarry Smith 
299d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCSetFromOptions_Redundant(PC pc, PetscOptionItems *PetscOptionsObject)
300d71ae5a4SJacob Faibussowitsch {
301a98ce0f4SHong Zhang   PC_Redundant *red = (PC_Redundant *)pc->data;
302a98ce0f4SHong Zhang 
3034b9ad928SBarry Smith   PetscFunctionBegin;
304d0609cedSBarry Smith   PetscOptionsHeadBegin(PetscOptionsObject, "Redundant options");
3059566063dSJacob Faibussowitsch   PetscCall(PetscOptionsInt("-pc_redundant_number", "Number of redundant pc", "PCRedundantSetNumber", red->nsubcomm, &red->nsubcomm, NULL));
306d0609cedSBarry Smith   PetscOptionsHeadEnd();
3073ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
3084b9ad928SBarry Smith }
3094b9ad928SBarry Smith 
310d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCRedundantSetNumber_Redundant(PC pc, PetscInt nreds)
311d71ae5a4SJacob Faibussowitsch {
31209a6bc64SHong Zhang   PC_Redundant *red = (PC_Redundant *)pc->data;
31309a6bc64SHong Zhang 
31409a6bc64SHong Zhang   PetscFunctionBegin;
31509a6bc64SHong Zhang   red->nsubcomm = nreds;
3163ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
31709a6bc64SHong Zhang }
31809a6bc64SHong Zhang 
31909a6bc64SHong Zhang /*@
32009a6bc64SHong Zhang    PCRedundantSetNumber - Sets the number of redundant preconditioner contexts.
32109a6bc64SHong Zhang 
322c3339decSBarry Smith    Logically Collective
32309a6bc64SHong Zhang 
32409a6bc64SHong Zhang    Input Parameters:
32509a6bc64SHong Zhang +  pc - the preconditioner context
3269b21d695SBarry Smith -  nredundant - number of redundant preconditioner contexts; for example if you are using 64 MPI processes and
3279b21d695SBarry Smith                               use an nredundant of 4 there will be 4 parallel solves each on 16 = 64/4 processes.
32809a6bc64SHong Zhang 
32909a6bc64SHong Zhang    Level: advanced
33009a6bc64SHong Zhang 
331f1580f4eSBarry Smith .seealso: `PCREDUNDANT`
33209a6bc64SHong Zhang @*/
333d71ae5a4SJacob Faibussowitsch PetscErrorCode PCRedundantSetNumber(PC pc, PetscInt nredundant)
334d71ae5a4SJacob Faibussowitsch {
33509a6bc64SHong Zhang   PetscFunctionBegin;
3360700a824SBarry Smith   PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
33763a3b9bcSJacob Faibussowitsch   PetscCheck(nredundant > 0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "num of redundant pc %" PetscInt_FMT " must be positive", nredundant);
338cac4c232SBarry Smith   PetscTryMethod(pc, "PCRedundantSetNumber_C", (PC, PetscInt), (pc, nredundant));
3393ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
34009a6bc64SHong Zhang }
34109a6bc64SHong Zhang 
342d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCRedundantSetScatter_Redundant(PC pc, VecScatter in, VecScatter out)
343d71ae5a4SJacob Faibussowitsch {
3444b9ad928SBarry Smith   PC_Redundant *red = (PC_Redundant *)pc->data;
3454b9ad928SBarry Smith 
3464b9ad928SBarry Smith   PetscFunctionBegin;
3479566063dSJacob Faibussowitsch   PetscCall(PetscObjectReference((PetscObject)in));
3489566063dSJacob Faibussowitsch   PetscCall(VecScatterDestroy(&red->scatterin));
3492fa5cd67SKarl Rupp 
350c3122656SLisandro Dalcin   red->scatterin = in;
3512fa5cd67SKarl Rupp 
3529566063dSJacob Faibussowitsch   PetscCall(PetscObjectReference((PetscObject)out));
3539566063dSJacob Faibussowitsch   PetscCall(VecScatterDestroy(&red->scatterout));
354c3122656SLisandro Dalcin   red->scatterout = out;
3553ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
3564b9ad928SBarry Smith }
3574b9ad928SBarry Smith 
3584b9ad928SBarry Smith /*@
3594b9ad928SBarry Smith    PCRedundantSetScatter - Sets the scatter used to copy values into the
3604b9ad928SBarry Smith      redundant local solve and the scatter to move them back into the global
3614b9ad928SBarry Smith      vector.
3624b9ad928SBarry Smith 
363c3339decSBarry Smith    Logically Collective
3644b9ad928SBarry Smith 
3654b9ad928SBarry Smith    Input Parameters:
3664b9ad928SBarry Smith +  pc - the preconditioner context
3674b9ad928SBarry Smith .  in - the scatter to move the values in
3684b9ad928SBarry Smith -  out - the scatter to move them out
3694b9ad928SBarry Smith 
3704b9ad928SBarry Smith    Level: advanced
3714b9ad928SBarry Smith 
372f1580f4eSBarry Smith .seealso: `PCREDUNDANT`
3734b9ad928SBarry Smith @*/
374d71ae5a4SJacob Faibussowitsch PetscErrorCode PCRedundantSetScatter(PC pc, VecScatter in, VecScatter out)
375d71ae5a4SJacob Faibussowitsch {
3764b9ad928SBarry Smith   PetscFunctionBegin;
3770700a824SBarry Smith   PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
37897929ea7SJunchao Zhang   PetscValidHeaderSpecific(in, PETSCSF_CLASSID, 2);
37997929ea7SJunchao Zhang   PetscValidHeaderSpecific(out, PETSCSF_CLASSID, 3);
380cac4c232SBarry Smith   PetscTryMethod(pc, "PCRedundantSetScatter_C", (PC, VecScatter, VecScatter), (pc, in, out));
3813ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
3824b9ad928SBarry Smith }
3834b9ad928SBarry Smith 
384d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCRedundantGetKSP_Redundant(PC pc, KSP *innerksp)
385d71ae5a4SJacob Faibussowitsch {
3864b9ad928SBarry Smith   PC_Redundant *red = (PC_Redundant *)pc->data;
38775024027SHong Zhang   MPI_Comm      comm, subcomm;
38875024027SHong Zhang   const char   *prefix;
38908cecb0aSPierre Jolivet   PetscBool     issbaij;
3904b9ad928SBarry Smith 
3914b9ad928SBarry Smith   PetscFunctionBegin;
39275024027SHong Zhang   if (!red->psubcomm) {
3939566063dSJacob Faibussowitsch     PetscCall(PCGetOptionsPrefix(pc, &prefix));
394e5acf8a4SHong Zhang 
3959566063dSJacob Faibussowitsch     PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
3969566063dSJacob Faibussowitsch     PetscCall(PetscSubcommCreate(comm, &red->psubcomm));
3979566063dSJacob Faibussowitsch     PetscCall(PetscSubcommSetNumber(red->psubcomm, red->nsubcomm));
3989566063dSJacob Faibussowitsch     PetscCall(PetscSubcommSetType(red->psubcomm, PETSC_SUBCOMM_CONTIGUOUS));
399e5acf8a4SHong Zhang 
4009566063dSJacob Faibussowitsch     PetscCall(PetscSubcommSetOptionsPrefix(red->psubcomm, prefix));
4019566063dSJacob Faibussowitsch     PetscCall(PetscSubcommSetFromOptions(red->psubcomm));
40275024027SHong Zhang 
40375024027SHong Zhang     /* create a new PC that processors in each subcomm have copy of */
40475024027SHong Zhang     subcomm = PetscSubcommChild(red->psubcomm);
40575024027SHong Zhang 
4069566063dSJacob Faibussowitsch     PetscCall(KSPCreate(subcomm, &red->ksp));
4079566063dSJacob Faibussowitsch     PetscCall(KSPSetErrorIfNotConverged(red->ksp, pc->erroriffailure));
4089566063dSJacob Faibussowitsch     PetscCall(PetscObjectIncrementTabLevel((PetscObject)red->ksp, (PetscObject)pc, 1));
4099566063dSJacob Faibussowitsch     PetscCall(KSPSetType(red->ksp, KSPPREONLY));
4109566063dSJacob Faibussowitsch     PetscCall(KSPGetPC(red->ksp, &red->pc));
4119566063dSJacob Faibussowitsch     PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATSEQSBAIJ, &issbaij));
41248a46eb9SPierre Jolivet     if (!issbaij) PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATMPISBAIJ, &issbaij));
41308cecb0aSPierre Jolivet     if (!issbaij) {
4149566063dSJacob Faibussowitsch       PetscCall(PCSetType(red->pc, PCLU));
41508cecb0aSPierre Jolivet     } else {
4169566063dSJacob Faibussowitsch       PetscCall(PCSetType(red->pc, PCCHOLESKY));
41708cecb0aSPierre Jolivet     }
418753b7fb9SBarry Smith     if (red->shifttypeset) {
4199566063dSJacob Faibussowitsch       PetscCall(PCFactorSetShiftType(red->pc, red->shifttype));
420753b7fb9SBarry Smith       red->shifttypeset = PETSC_FALSE;
421753b7fb9SBarry Smith     }
4229566063dSJacob Faibussowitsch     PetscCall(KSPSetOptionsPrefix(red->ksp, prefix));
4239566063dSJacob Faibussowitsch     PetscCall(KSPAppendOptionsPrefix(red->ksp, "redundant_"));
42475024027SHong Zhang   }
42583ab6a24SBarry Smith   *innerksp = red->ksp;
4263ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
4274b9ad928SBarry Smith }
4284b9ad928SBarry Smith 
4294b9ad928SBarry Smith /*@
430f1580f4eSBarry Smith    PCRedundantGetKSP - Gets the less parallel `KSP` created by the redundant `PC`.
4314b9ad928SBarry Smith 
4324b9ad928SBarry Smith    Not Collective
4334b9ad928SBarry Smith 
4344b9ad928SBarry Smith    Input Parameter:
4354b9ad928SBarry Smith .  pc - the preconditioner context
4364b9ad928SBarry Smith 
4374b9ad928SBarry Smith    Output Parameter:
438f1580f4eSBarry Smith .  innerksp - the `KSP` on the smaller set of processes
4394b9ad928SBarry Smith 
4404b9ad928SBarry Smith    Level: advanced
4414b9ad928SBarry Smith 
442f1580f4eSBarry Smith .seealso: `PCREDUNDANT`
4434b9ad928SBarry Smith @*/
444d71ae5a4SJacob Faibussowitsch PetscErrorCode PCRedundantGetKSP(PC pc, KSP *innerksp)
445d71ae5a4SJacob Faibussowitsch {
4464b9ad928SBarry Smith   PetscFunctionBegin;
4470700a824SBarry Smith   PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
44883ab6a24SBarry Smith   PetscValidPointer(innerksp, 2);
449cac4c232SBarry Smith   PetscUseMethod(pc, "PCRedundantGetKSP_C", (PC, KSP *), (pc, innerksp));
4503ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
4514b9ad928SBarry Smith }
4524b9ad928SBarry Smith 
453d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCRedundantGetOperators_Redundant(PC pc, Mat *mat, Mat *pmat)
454d71ae5a4SJacob Faibussowitsch {
4554b9ad928SBarry Smith   PC_Redundant *red = (PC_Redundant *)pc->data;
4564b9ad928SBarry Smith 
4574b9ad928SBarry Smith   PetscFunctionBegin;
458b3804887SHong Zhang   if (mat) *mat = red->pmats;
459b3804887SHong Zhang   if (pmat) *pmat = red->pmats;
4603ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
4614b9ad928SBarry Smith }
4624b9ad928SBarry Smith 
4634b9ad928SBarry Smith /*@
4644b9ad928SBarry Smith    PCRedundantGetOperators - gets the sequential matrix and preconditioner matrix
4654b9ad928SBarry Smith 
4664b9ad928SBarry Smith    Not Collective
4674b9ad928SBarry Smith 
4684b9ad928SBarry Smith    Input Parameter:
4694b9ad928SBarry Smith .  pc - the preconditioner context
4704b9ad928SBarry Smith 
4714b9ad928SBarry Smith    Output Parameters:
4724b9ad928SBarry Smith +  mat - the matrix
4734b9ad928SBarry Smith -  pmat - the (possibly different) preconditioner matrix
4744b9ad928SBarry Smith 
4754b9ad928SBarry Smith    Level: advanced
4764b9ad928SBarry Smith 
477f1580f4eSBarry Smith .seealso: `PCREDUNDANT`
4784b9ad928SBarry Smith @*/
479d71ae5a4SJacob Faibussowitsch PetscErrorCode PCRedundantGetOperators(PC pc, Mat *mat, Mat *pmat)
480d71ae5a4SJacob Faibussowitsch {
4814b9ad928SBarry Smith   PetscFunctionBegin;
4820700a824SBarry Smith   PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
4834482741eSBarry Smith   if (mat) PetscValidPointer(mat, 2);
4844482741eSBarry Smith   if (pmat) PetscValidPointer(pmat, 3);
485cac4c232SBarry Smith   PetscUseMethod(pc, "PCRedundantGetOperators_C", (PC, Mat *, Mat *), (pc, mat, pmat));
4863ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
4874b9ad928SBarry Smith }
4884b9ad928SBarry Smith 
48937a17b4dSBarry Smith /*MC
490f1580f4eSBarry Smith      PCREDUNDANT - Runs a `KSP` solver with preconditioner for the entire problem on subgroups of processors
49137a17b4dSBarry Smith 
492f1580f4eSBarry Smith   Options Database Key:
4939b21d695SBarry Smith .  -pc_redundant_number <n> - number of redundant solves, for example if you are using 64 MPI processes and
4949b21d695SBarry Smith                               use an n of 4 there will be 4 parallel solves each on 16 = 64/4 processes.
49509391456SBarry Smith 
49637a17b4dSBarry Smith    Level: intermediate
49737a17b4dSBarry Smith 
49895452b02SPatrick Sanan    Notes:
499f1580f4eSBarry Smith    Options for the redundant preconditioners can be set using the options database prefix -redundant_
50083ab6a24SBarry Smith 
501f1580f4eSBarry Smith    The default `KSP` is preonly and the default `PC` is `PCLU` or `PCCHOLESKY` if Pmat is of type `MATSBAIJ`.
502753b7fb9SBarry Smith 
503f1580f4eSBarry Smith    `PCFactorSetShiftType()` applied to this `PC` will convey they shift type into the inner `PC` if it is factorization based.
504f1580f4eSBarry Smith 
505f1580f4eSBarry Smith    Developer Note:
506f1580f4eSBarry Smith    `PCSetInitialGuessNonzero()` is not used by this class but likely should be.
5079cfaa89bSBarry Smith 
508db781477SPatrick Sanan .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PCRedundantSetScatter()`,
509f1580f4eSBarry Smith           `PCRedundantGetKSP()`, `PCRedundantGetOperators()`, `PCRedundantSetNumber()`, `PCREDISTRIBUTE`
51037a17b4dSBarry Smith M*/
51137a17b4dSBarry Smith 
512d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode PCCreate_Redundant(PC pc)
513d71ae5a4SJacob Faibussowitsch {
5144b9ad928SBarry Smith   PC_Redundant *red;
51569db28dcSHong Zhang   PetscMPIInt   size;
5163f457be1SHong Zhang 
5174b9ad928SBarry Smith   PetscFunctionBegin;
5184dfa11a4SJacob Faibussowitsch   PetscCall(PetscNew(&red));
5199566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
5202fa5cd67SKarl Rupp 
52169db28dcSHong Zhang   red->nsubcomm       = size;
5224b9ad928SBarry Smith   red->useparallelmat = PETSC_TRUE;
5231fbd8f88SHong Zhang   pc->data            = (void *)red;
5244b9ad928SBarry Smith 
5254b9ad928SBarry Smith   pc->ops->apply          = PCApply_Redundant;
526d88bfacbSStefano Zampini   pc->ops->applytranspose = PCApplyTranspose_Redundant;
5274b9ad928SBarry Smith   pc->ops->setup          = PCSetUp_Redundant;
5284b9ad928SBarry Smith   pc->ops->destroy        = PCDestroy_Redundant;
5291ea5a559SBarry Smith   pc->ops->reset          = PCReset_Redundant;
5304b9ad928SBarry Smith   pc->ops->setfromoptions = PCSetFromOptions_Redundant;
5314b9ad928SBarry Smith   pc->ops->view           = PCView_Redundant;
5322fa5cd67SKarl Rupp 
5339566063dSJacob Faibussowitsch   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantSetScatter_C", PCRedundantSetScatter_Redundant));
5349566063dSJacob Faibussowitsch   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantSetNumber_C", PCRedundantSetNumber_Redundant));
5359566063dSJacob Faibussowitsch   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantGetKSP_C", PCRedundantGetKSP_Redundant));
5369566063dSJacob Faibussowitsch   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantGetOperators_C", PCRedundantGetOperators_Redundant));
5379566063dSJacob Faibussowitsch   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFactorSetShiftType_C", PCFactorSetShiftType_Redundant));
5383ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
5394b9ad928SBarry Smith }
540