14b9ad928SBarry Smith /*
23f457be1SHong Zhang This file defines a "solve the problem redundantly on each subgroup of processor" preconditioner.
34b9ad928SBarry Smith */
4af0996ceSBarry Smith #include <petsc/private/pcimpl.h>
5c6db04a5SJed Brown #include <petscksp.h> /*I "petscksp.h" I*/
64b9ad928SBarry Smith
74b9ad928SBarry Smith typedef struct {
83e065800SHong Zhang KSP ksp;
94b9ad928SBarry Smith PC pc; /* actual preconditioner used on each processor */
10ce94432eSBarry Smith Vec xsub, ysub; /* vectors of a subcommunicator to hold parallel vectors of PetscObjectComm((PetscObject)pc) */
113f457be1SHong Zhang Vec xdup, ydup; /* parallel vector that congregates xsub or ysub facilitating vector scattering */
127addb90fSBarry Smith Mat pmats; /* matrices belong to a subcommunicator */
133f457be1SHong Zhang VecScatter scatterin, scatterout; /* scatter used to move all values to each processor group (subcommunicator) */
14ace3abfcSBarry Smith PetscBool useparallelmat;
15c540e29cSHong Zhang PetscSubcomm psubcomm;
161fbd8f88SHong Zhang PetscInt nsubcomm; /* num of data structure PetscSubcomm */
17753b7fb9SBarry Smith PetscBool shifttypeset;
18753b7fb9SBarry Smith MatFactorShiftType shifttype;
194b9ad928SBarry Smith } PC_Redundant;
204b9ad928SBarry Smith
PCFactorSetShiftType_Redundant(PC pc,MatFactorShiftType shifttype)2166976f2fSJacob Faibussowitsch static PetscErrorCode PCFactorSetShiftType_Redundant(PC pc, MatFactorShiftType shifttype)
22d71ae5a4SJacob Faibussowitsch {
23753b7fb9SBarry Smith PC_Redundant *red = (PC_Redundant *)pc->data;
24753b7fb9SBarry Smith
25753b7fb9SBarry Smith PetscFunctionBegin;
26753b7fb9SBarry Smith if (red->ksp) {
27753b7fb9SBarry Smith PC pc;
289566063dSJacob Faibussowitsch PetscCall(KSPGetPC(red->ksp, &pc));
299566063dSJacob Faibussowitsch PetscCall(PCFactorSetShiftType(pc, shifttype));
30753b7fb9SBarry Smith } else {
31753b7fb9SBarry Smith red->shifttypeset = PETSC_TRUE;
32753b7fb9SBarry Smith red->shifttype = shifttype;
33753b7fb9SBarry Smith }
343ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
35753b7fb9SBarry Smith }
36753b7fb9SBarry Smith
PCView_Redundant(PC pc,PetscViewer viewer)37d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCView_Redundant(PC pc, PetscViewer viewer)
38d71ae5a4SJacob Faibussowitsch {
394b9ad928SBarry Smith PC_Redundant *red = (PC_Redundant *)pc->data;
40*9f196a02SMartin Diehl PetscBool isascii, isstring;
4103ccd0b4SBarry Smith PetscViewer subviewer;
424b9ad928SBarry Smith
434b9ad928SBarry Smith PetscFunctionBegin;
44*9f196a02SMartin Diehl PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
459566063dSJacob Faibussowitsch PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
46*9f196a02SMartin Diehl if (isascii) {
4703ccd0b4SBarry Smith if (!red->psubcomm) {
489566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " Not yet setup\n"));
4903ccd0b4SBarry Smith } else {
5063a3b9bcSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " First (color=0) of %" PetscInt_FMT " PCs follows\n", red->nsubcomm));
519566063dSJacob Faibussowitsch PetscCall(PetscViewerGetSubViewer(viewer, ((PetscObject)red->pc)->comm, &subviewer));
52f5dd71faSBarry Smith if (!red->psubcomm->color) { /* only view first redundant pc */
539566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPushTab(subviewer));
549566063dSJacob Faibussowitsch PetscCall(KSPView(red->ksp, subviewer));
559566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPopTab(subviewer));
564b9ad928SBarry Smith }
579566063dSJacob Faibussowitsch PetscCall(PetscViewerRestoreSubViewer(viewer, ((PetscObject)red->pc)->comm, &subviewer));
584b9ad928SBarry Smith }
5903ccd0b4SBarry Smith } else if (isstring) {
609566063dSJacob Faibussowitsch PetscCall(PetscViewerStringSPrintf(viewer, " Redundant solver preconditioner"));
614b9ad928SBarry Smith }
623ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
634b9ad928SBarry Smith }
644b9ad928SBarry Smith
6519b3b6edSHong Zhang #include <../src/mat/impls/aij/mpi/mpiaij.h>
PCSetUp_Redundant(PC pc)66d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCSetUp_Redundant(PC pc)
67d71ae5a4SJacob Faibussowitsch {
684b9ad928SBarry Smith PC_Redundant *red = (PC_Redundant *)pc->data;
691b81debcSHong Zhang PetscInt mstart, mend, mlocal, M;
7013f74950SBarry Smith PetscMPIInt size;
71ce94432eSBarry Smith MPI_Comm comm, subcomm;
72ddc54837SHong Zhang Vec x;
733f457be1SHong Zhang
744b9ad928SBarry Smith PetscFunctionBegin;
759566063dSJacob Faibussowitsch PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
76ddc54837SHong Zhang
77ddc54837SHong Zhang /* if pmatrix set by user is sequential then we do not need to gather the parallel matrix */
789566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_size(comm, &size));
79ddc54837SHong Zhang if (size == 1) red->useparallelmat = PETSC_FALSE;
801fbd8f88SHong Zhang
814b9ad928SBarry Smith if (!pc->setupcalled) {
821b81debcSHong Zhang PetscInt mloc_sub;
8375024027SHong Zhang if (!red->psubcomm) { /* create red->psubcomm, new ksp and pc over subcomm */
8475024027SHong Zhang KSP ksp;
859566063dSJacob Faibussowitsch PetscCall(PCRedundantGetKSP(pc, &ksp));
861b81debcSHong Zhang }
8775024027SHong Zhang subcomm = PetscSubcommChild(red->psubcomm);
881fbd8f88SHong Zhang
891b81debcSHong Zhang if (red->useparallelmat) {
90268865b1SPierre Jolivet /* grab the parallel matrix and put it into the processes of a subcommunicator */
919566063dSJacob Faibussowitsch PetscCall(MatCreateRedundantMatrix(pc->pmat, red->psubcomm->n, subcomm, MAT_INITIAL_MATRIX, &red->pmats));
92b85f2e9bSHong Zhang
939566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_size(subcomm, &size));
94b85f2e9bSHong Zhang if (size > 1) {
9508cecb0aSPierre Jolivet PetscBool foundpack, issbaij;
969566063dSJacob Faibussowitsch PetscCall(PetscObjectTypeCompare((PetscObject)red->pmats, MATMPISBAIJ, &issbaij));
9708cecb0aSPierre Jolivet if (!issbaij) {
989566063dSJacob Faibussowitsch PetscCall(MatGetFactorAvailable(red->pmats, NULL, MAT_FACTOR_LU, &foundpack));
9908cecb0aSPierre Jolivet } else {
1009566063dSJacob Faibussowitsch PetscCall(MatGetFactorAvailable(red->pmats, NULL, MAT_FACTOR_CHOLESKY, &foundpack));
10108cecb0aSPierre Jolivet }
102b85f2e9bSHong Zhang if (!foundpack) { /* reset default ksp and pc */
1039566063dSJacob Faibussowitsch PetscCall(KSPSetType(red->ksp, KSPGMRES));
1049566063dSJacob Faibussowitsch PetscCall(PCSetType(red->pc, PCBJACOBI));
105c1619fb6SBarry Smith } else {
1069566063dSJacob Faibussowitsch PetscCall(PCFactorSetMatSolverType(red->pc, NULL));
107b85f2e9bSHong Zhang }
108b85f2e9bSHong Zhang }
109b85f2e9bSHong Zhang
1109566063dSJacob Faibussowitsch PetscCall(KSPSetOperators(red->ksp, red->pmats, red->pmats));
1114b9ad928SBarry Smith
1121b81debcSHong Zhang /* get working vectors xsub and ysub */
1139566063dSJacob Faibussowitsch PetscCall(MatCreateVecs(red->pmats, &red->xsub, &red->ysub));
1142fa5cd67SKarl Rupp
1158b96b0d2SHong Zhang /* create working vectors xdup and ydup.
11646091a0eSPierre Jolivet xdup concatenates all xsub's contiguously to form a mpi vector over dupcomm (see PetscSubcommCreate_interlaced())
1178b96b0d2SHong Zhang ydup concatenates all ysub and has empty local arrays because ysub's arrays will be place into it.
118ce94432eSBarry Smith Note: we use communicator dupcomm, not PetscObjectComm((PetscObject)pc)! */
1199566063dSJacob Faibussowitsch PetscCall(MatGetLocalSize(red->pmats, &mloc_sub, NULL));
1209566063dSJacob Faibussowitsch PetscCall(VecCreateMPI(PetscSubcommContiguousParent(red->psubcomm), mloc_sub, PETSC_DECIDE, &red->xdup));
1219566063dSJacob Faibussowitsch PetscCall(VecCreateMPIWithArray(PetscSubcommContiguousParent(red->psubcomm), 1, mloc_sub, PETSC_DECIDE, NULL, &red->ydup));
1223f457be1SHong Zhang
123f68be91cSHong Zhang /* create vecscatters */
124f68be91cSHong Zhang if (!red->scatterin) { /* efficiency of scatterin is independent from psubcomm_type! */
1253f457be1SHong Zhang IS is1, is2;
1263f457be1SHong Zhang PetscInt *idx1, *idx2, i, j, k;
12745fc02eaSBarry Smith
1289566063dSJacob Faibussowitsch PetscCall(MatCreateVecs(pc->pmat, &x, NULL));
1299566063dSJacob Faibussowitsch PetscCall(VecGetSize(x, &M));
1309566063dSJacob Faibussowitsch PetscCall(VecGetOwnershipRange(x, &mstart, &mend));
1311b81debcSHong Zhang mlocal = mend - mstart;
1329566063dSJacob Faibussowitsch PetscCall(PetscMalloc2(red->psubcomm->n * mlocal, &idx1, red->psubcomm->n * mlocal, &idx2));
1333f457be1SHong Zhang j = 0;
1341fbd8f88SHong Zhang for (k = 0; k < red->psubcomm->n; k++) {
1353f457be1SHong Zhang for (i = mstart; i < mend; i++) {
1363f457be1SHong Zhang idx1[j] = i;
137ddc54837SHong Zhang idx2[j++] = i + M * k;
1383f457be1SHong Zhang }
1393f457be1SHong Zhang }
1409566063dSJacob Faibussowitsch PetscCall(ISCreateGeneral(comm, red->psubcomm->n * mlocal, idx1, PETSC_COPY_VALUES, &is1));
1419566063dSJacob Faibussowitsch PetscCall(ISCreateGeneral(comm, red->psubcomm->n * mlocal, idx2, PETSC_COPY_VALUES, &is2));
1429566063dSJacob Faibussowitsch PetscCall(VecScatterCreate(x, is1, red->xdup, is2, &red->scatterin));
1439566063dSJacob Faibussowitsch PetscCall(ISDestroy(&is1));
1449566063dSJacob Faibussowitsch PetscCall(ISDestroy(&is2));
1453f457be1SHong Zhang
1466909748bSHong Zhang /* Impl below is good for PETSC_SUBCOMM_INTERLACED (no inter-process communication) and PETSC_SUBCOMM_CONTIGUOUS (communication within subcomm) */
1479566063dSJacob Faibussowitsch PetscCall(ISCreateStride(comm, mlocal, mstart + red->psubcomm->color * M, 1, &is1));
1489566063dSJacob Faibussowitsch PetscCall(ISCreateStride(comm, mlocal, mstart, 1, &is2));
1499566063dSJacob Faibussowitsch PetscCall(VecScatterCreate(red->xdup, is1, x, is2, &red->scatterout));
1509566063dSJacob Faibussowitsch PetscCall(ISDestroy(&is1));
1519566063dSJacob Faibussowitsch PetscCall(ISDestroy(&is2));
1529566063dSJacob Faibussowitsch PetscCall(PetscFree2(idx1, idx2));
1539566063dSJacob Faibussowitsch PetscCall(VecDestroy(&x));
1541b81debcSHong Zhang }
155ab661555SHong Zhang } else { /* !red->useparallelmat */
1569566063dSJacob Faibussowitsch PetscCall(KSPSetOperators(red->ksp, pc->mat, pc->pmat));
1571b81debcSHong Zhang }
158ab661555SHong Zhang } else { /* pc->setupcalled */
1594b9ad928SBarry Smith if (red->useparallelmat) {
160ab661555SHong Zhang MatReuse reuse;
161268865b1SPierre Jolivet /* grab the parallel matrix and put it into the processes of a subcommunicator */
162ab661555SHong Zhang if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1634b9ad928SBarry Smith /* destroy old matrices */
1649566063dSJacob Faibussowitsch PetscCall(MatDestroy(&red->pmats));
165ab661555SHong Zhang reuse = MAT_INITIAL_MATRIX;
1664b9ad928SBarry Smith } else {
167ab661555SHong Zhang reuse = MAT_REUSE_MATRIX;
168ab661555SHong Zhang }
1699566063dSJacob Faibussowitsch PetscCall(MatCreateRedundantMatrix(pc->pmat, red->psubcomm->n, PetscSubcommChild(red->psubcomm), reuse, &red->pmats));
1709566063dSJacob Faibussowitsch PetscCall(KSPSetOperators(red->ksp, red->pmats, red->pmats));
171ab661555SHong Zhang } else { /* !red->useparallelmat */
1729566063dSJacob Faibussowitsch PetscCall(KSPSetOperators(red->ksp, pc->mat, pc->pmat));
1734b9ad928SBarry Smith }
174ab661555SHong Zhang }
1751b81debcSHong Zhang
1761baa6e33SBarry Smith if (pc->setfromoptionscalled) PetscCall(KSPSetFromOptions(red->ksp));
1773ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
1784b9ad928SBarry Smith }
1794b9ad928SBarry Smith
PCApply_Redundant(PC pc,Vec x,Vec y)180d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCApply_Redundant(PC pc, Vec x, Vec y)
181d71ae5a4SJacob Faibussowitsch {
1824b9ad928SBarry Smith PC_Redundant *red = (PC_Redundant *)pc->data;
1833f457be1SHong Zhang PetscScalar *array;
1844b9ad928SBarry Smith
1854b9ad928SBarry Smith PetscFunctionBegin;
186ddc54837SHong Zhang if (!red->useparallelmat) {
1879566063dSJacob Faibussowitsch PetscCall(KSPSolve(red->ksp, x, y));
1889566063dSJacob Faibussowitsch PetscCall(KSPCheckSolve(red->ksp, pc, y));
1893ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
190ddc54837SHong Zhang }
191ddc54837SHong Zhang
1923f457be1SHong Zhang /* scatter x to xdup */
1939566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(red->scatterin, x, red->xdup, INSERT_VALUES, SCATTER_FORWARD));
1949566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(red->scatterin, x, red->xdup, INSERT_VALUES, SCATTER_FORWARD));
1953f457be1SHong Zhang
1963f457be1SHong Zhang /* place xdup's local array into xsub */
1979566063dSJacob Faibussowitsch PetscCall(VecGetArray(red->xdup, &array));
1989566063dSJacob Faibussowitsch PetscCall(VecPlaceArray(red->xsub, (const PetscScalar *)array));
1994b9ad928SBarry Smith
2004b9ad928SBarry Smith /* apply preconditioner on each processor */
2019566063dSJacob Faibussowitsch PetscCall(KSPSolve(red->ksp, red->xsub, red->ysub));
2029566063dSJacob Faibussowitsch PetscCall(KSPCheckSolve(red->ksp, pc, red->ysub));
2039566063dSJacob Faibussowitsch PetscCall(VecResetArray(red->xsub));
2049566063dSJacob Faibussowitsch PetscCall(VecRestoreArray(red->xdup, &array));
2054b9ad928SBarry Smith
2063f457be1SHong Zhang /* place ysub's local array into ydup */
2079566063dSJacob Faibussowitsch PetscCall(VecGetArray(red->ysub, &array));
2089566063dSJacob Faibussowitsch PetscCall(VecPlaceArray(red->ydup, (const PetscScalar *)array));
2093f457be1SHong Zhang
2103f457be1SHong Zhang /* scatter ydup to y */
2119566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(red->scatterout, red->ydup, y, INSERT_VALUES, SCATTER_FORWARD));
2129566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(red->scatterout, red->ydup, y, INSERT_VALUES, SCATTER_FORWARD));
2139566063dSJacob Faibussowitsch PetscCall(VecResetArray(red->ydup));
2149566063dSJacob Faibussowitsch PetscCall(VecRestoreArray(red->ysub, &array));
2153ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
2164b9ad928SBarry Smith }
2174b9ad928SBarry Smith
PCApplyTranspose_Redundant(PC pc,Vec x,Vec y)218d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCApplyTranspose_Redundant(PC pc, Vec x, Vec y)
219d71ae5a4SJacob Faibussowitsch {
220d88bfacbSStefano Zampini PC_Redundant *red = (PC_Redundant *)pc->data;
221d88bfacbSStefano Zampini PetscScalar *array;
222d88bfacbSStefano Zampini
223d88bfacbSStefano Zampini PetscFunctionBegin;
224d88bfacbSStefano Zampini if (!red->useparallelmat) {
2259566063dSJacob Faibussowitsch PetscCall(KSPSolveTranspose(red->ksp, x, y));
2269566063dSJacob Faibussowitsch PetscCall(KSPCheckSolve(red->ksp, pc, y));
2273ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
228d88bfacbSStefano Zampini }
229d88bfacbSStefano Zampini
230d88bfacbSStefano Zampini /* scatter x to xdup */
2319566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(red->scatterin, x, red->xdup, INSERT_VALUES, SCATTER_FORWARD));
2329566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(red->scatterin, x, red->xdup, INSERT_VALUES, SCATTER_FORWARD));
233d88bfacbSStefano Zampini
234d88bfacbSStefano Zampini /* place xdup's local array into xsub */
2359566063dSJacob Faibussowitsch PetscCall(VecGetArray(red->xdup, &array));
2369566063dSJacob Faibussowitsch PetscCall(VecPlaceArray(red->xsub, (const PetscScalar *)array));
237d88bfacbSStefano Zampini
238d88bfacbSStefano Zampini /* apply preconditioner on each processor */
2399566063dSJacob Faibussowitsch PetscCall(KSPSolveTranspose(red->ksp, red->xsub, red->ysub));
2409566063dSJacob Faibussowitsch PetscCall(KSPCheckSolve(red->ksp, pc, red->ysub));
2419566063dSJacob Faibussowitsch PetscCall(VecResetArray(red->xsub));
2429566063dSJacob Faibussowitsch PetscCall(VecRestoreArray(red->xdup, &array));
243d88bfacbSStefano Zampini
244d88bfacbSStefano Zampini /* place ysub's local array into ydup */
2459566063dSJacob Faibussowitsch PetscCall(VecGetArray(red->ysub, &array));
2469566063dSJacob Faibussowitsch PetscCall(VecPlaceArray(red->ydup, (const PetscScalar *)array));
247d88bfacbSStefano Zampini
248d88bfacbSStefano Zampini /* scatter ydup to y */
2499566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(red->scatterout, red->ydup, y, INSERT_VALUES, SCATTER_FORWARD));
2509566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(red->scatterout, red->ydup, y, INSERT_VALUES, SCATTER_FORWARD));
2519566063dSJacob Faibussowitsch PetscCall(VecResetArray(red->ydup));
2529566063dSJacob Faibussowitsch PetscCall(VecRestoreArray(red->ysub, &array));
2533ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
254d88bfacbSStefano Zampini }
255d88bfacbSStefano Zampini
PCReset_Redundant(PC pc)256d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCReset_Redundant(PC pc)
257d71ae5a4SJacob Faibussowitsch {
2584b9ad928SBarry Smith PC_Redundant *red = (PC_Redundant *)pc->data;
2594b9ad928SBarry Smith
2604b9ad928SBarry Smith PetscFunctionBegin;
2611b81debcSHong Zhang if (red->useparallelmat) {
2629566063dSJacob Faibussowitsch PetscCall(VecScatterDestroy(&red->scatterin));
2639566063dSJacob Faibussowitsch PetscCall(VecScatterDestroy(&red->scatterout));
2649566063dSJacob Faibussowitsch PetscCall(VecDestroy(&red->ysub));
2659566063dSJacob Faibussowitsch PetscCall(VecDestroy(&red->xsub));
2669566063dSJacob Faibussowitsch PetscCall(VecDestroy(&red->xdup));
2679566063dSJacob Faibussowitsch PetscCall(VecDestroy(&red->ydup));
2681b81debcSHong Zhang }
2699566063dSJacob Faibussowitsch PetscCall(MatDestroy(&red->pmats));
2709566063dSJacob Faibussowitsch PetscCall(KSPReset(red->ksp));
2713ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
2721ea5a559SBarry Smith }
2731ea5a559SBarry Smith
PCDestroy_Redundant(PC pc)274d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCDestroy_Redundant(PC pc)
275d71ae5a4SJacob Faibussowitsch {
2761ea5a559SBarry Smith PC_Redundant *red = (PC_Redundant *)pc->data;
2771ea5a559SBarry Smith
2781ea5a559SBarry Smith PetscFunctionBegin;
2799566063dSJacob Faibussowitsch PetscCall(PCReset_Redundant(pc));
2809566063dSJacob Faibussowitsch PetscCall(KSPDestroy(&red->ksp));
2819566063dSJacob Faibussowitsch PetscCall(PetscSubcommDestroy(&red->psubcomm));
2822e956fe4SStefano Zampini PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantSetScatter_C", NULL));
2832e956fe4SStefano Zampini PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantSetNumber_C", NULL));
2842e956fe4SStefano Zampini PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantGetKSP_C", NULL));
2852e956fe4SStefano Zampini PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantGetOperators_C", NULL));
2862e956fe4SStefano Zampini PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFactorSetShiftType_C", NULL));
2879566063dSJacob Faibussowitsch PetscCall(PetscFree(pc->data));
2883ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
2894b9ad928SBarry Smith }
2904b9ad928SBarry Smith
PCSetFromOptions_Redundant(PC pc,PetscOptionItems PetscOptionsObject)291ce78bad3SBarry Smith static PetscErrorCode PCSetFromOptions_Redundant(PC pc, PetscOptionItems PetscOptionsObject)
292d71ae5a4SJacob Faibussowitsch {
293a98ce0f4SHong Zhang PC_Redundant *red = (PC_Redundant *)pc->data;
294a98ce0f4SHong Zhang
2954b9ad928SBarry Smith PetscFunctionBegin;
296d0609cedSBarry Smith PetscOptionsHeadBegin(PetscOptionsObject, "Redundant options");
2979566063dSJacob Faibussowitsch PetscCall(PetscOptionsInt("-pc_redundant_number", "Number of redundant pc", "PCRedundantSetNumber", red->nsubcomm, &red->nsubcomm, NULL));
298d0609cedSBarry Smith PetscOptionsHeadEnd();
2993ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
3004b9ad928SBarry Smith }
3014b9ad928SBarry Smith
PCRedundantSetNumber_Redundant(PC pc,PetscInt nreds)302d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCRedundantSetNumber_Redundant(PC pc, PetscInt nreds)
303d71ae5a4SJacob Faibussowitsch {
30409a6bc64SHong Zhang PC_Redundant *red = (PC_Redundant *)pc->data;
30509a6bc64SHong Zhang
30609a6bc64SHong Zhang PetscFunctionBegin;
30709a6bc64SHong Zhang red->nsubcomm = nreds;
3083ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
30909a6bc64SHong Zhang }
31009a6bc64SHong Zhang
31109a6bc64SHong Zhang /*@
31209a6bc64SHong Zhang PCRedundantSetNumber - Sets the number of redundant preconditioner contexts.
31309a6bc64SHong Zhang
314c3339decSBarry Smith Logically Collective
31509a6bc64SHong Zhang
31609a6bc64SHong Zhang Input Parameters:
31709a6bc64SHong Zhang + pc - the preconditioner context
3189b21d695SBarry Smith - nredundant - number of redundant preconditioner contexts; for example if you are using 64 MPI processes and
3199b21d695SBarry Smith use an nredundant of 4 there will be 4 parallel solves each on 16 = 64/4 processes.
32009a6bc64SHong Zhang
32109a6bc64SHong Zhang Level: advanced
32209a6bc64SHong Zhang
323562efe2eSBarry Smith .seealso: [](ch_ksp), `PCREDUNDANT`
32409a6bc64SHong Zhang @*/
PCRedundantSetNumber(PC pc,PetscInt nredundant)325d71ae5a4SJacob Faibussowitsch PetscErrorCode PCRedundantSetNumber(PC pc, PetscInt nredundant)
326d71ae5a4SJacob Faibussowitsch {
32709a6bc64SHong Zhang PetscFunctionBegin;
3280700a824SBarry Smith PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
32963a3b9bcSJacob Faibussowitsch PetscCheck(nredundant > 0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "num of redundant pc %" PetscInt_FMT " must be positive", nredundant);
330cac4c232SBarry Smith PetscTryMethod(pc, "PCRedundantSetNumber_C", (PC, PetscInt), (pc, nredundant));
3313ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
33209a6bc64SHong Zhang }
33309a6bc64SHong Zhang
PCRedundantSetScatter_Redundant(PC pc,VecScatter in,VecScatter out)334d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCRedundantSetScatter_Redundant(PC pc, VecScatter in, VecScatter out)
335d71ae5a4SJacob Faibussowitsch {
3364b9ad928SBarry Smith PC_Redundant *red = (PC_Redundant *)pc->data;
3374b9ad928SBarry Smith
3384b9ad928SBarry Smith PetscFunctionBegin;
3399566063dSJacob Faibussowitsch PetscCall(PetscObjectReference((PetscObject)in));
3409566063dSJacob Faibussowitsch PetscCall(VecScatterDestroy(&red->scatterin));
3412fa5cd67SKarl Rupp
342c3122656SLisandro Dalcin red->scatterin = in;
3432fa5cd67SKarl Rupp
3449566063dSJacob Faibussowitsch PetscCall(PetscObjectReference((PetscObject)out));
3459566063dSJacob Faibussowitsch PetscCall(VecScatterDestroy(&red->scatterout));
346c3122656SLisandro Dalcin red->scatterout = out;
3473ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
3484b9ad928SBarry Smith }
3494b9ad928SBarry Smith
3504b9ad928SBarry Smith /*@
3514b9ad928SBarry Smith PCRedundantSetScatter - Sets the scatter used to copy values into the
3524b9ad928SBarry Smith redundant local solve and the scatter to move them back into the global
3534b9ad928SBarry Smith vector.
3544b9ad928SBarry Smith
355c3339decSBarry Smith Logically Collective
3564b9ad928SBarry Smith
3574b9ad928SBarry Smith Input Parameters:
3584b9ad928SBarry Smith + pc - the preconditioner context
3594b9ad928SBarry Smith . in - the scatter to move the values in
3604b9ad928SBarry Smith - out - the scatter to move them out
3614b9ad928SBarry Smith
3624b9ad928SBarry Smith Level: advanced
3634b9ad928SBarry Smith
364562efe2eSBarry Smith .seealso: [](ch_ksp), `PCREDUNDANT`
3654b9ad928SBarry Smith @*/
PCRedundantSetScatter(PC pc,VecScatter in,VecScatter out)366d71ae5a4SJacob Faibussowitsch PetscErrorCode PCRedundantSetScatter(PC pc, VecScatter in, VecScatter out)
367d71ae5a4SJacob Faibussowitsch {
3684b9ad928SBarry Smith PetscFunctionBegin;
3690700a824SBarry Smith PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
37097929ea7SJunchao Zhang PetscValidHeaderSpecific(in, PETSCSF_CLASSID, 2);
37197929ea7SJunchao Zhang PetscValidHeaderSpecific(out, PETSCSF_CLASSID, 3);
372cac4c232SBarry Smith PetscTryMethod(pc, "PCRedundantSetScatter_C", (PC, VecScatter, VecScatter), (pc, in, out));
3733ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
3744b9ad928SBarry Smith }
3754b9ad928SBarry Smith
PCRedundantGetKSP_Redundant(PC pc,KSP * innerksp)376d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCRedundantGetKSP_Redundant(PC pc, KSP *innerksp)
377d71ae5a4SJacob Faibussowitsch {
3784b9ad928SBarry Smith PC_Redundant *red = (PC_Redundant *)pc->data;
37975024027SHong Zhang MPI_Comm comm, subcomm;
38075024027SHong Zhang const char *prefix;
38108cecb0aSPierre Jolivet PetscBool issbaij;
3824b9ad928SBarry Smith
3834b9ad928SBarry Smith PetscFunctionBegin;
38475024027SHong Zhang if (!red->psubcomm) {
3859566063dSJacob Faibussowitsch PetscCall(PCGetOptionsPrefix(pc, &prefix));
386e5acf8a4SHong Zhang
3879566063dSJacob Faibussowitsch PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
3889566063dSJacob Faibussowitsch PetscCall(PetscSubcommCreate(comm, &red->psubcomm));
3899566063dSJacob Faibussowitsch PetscCall(PetscSubcommSetNumber(red->psubcomm, red->nsubcomm));
3909566063dSJacob Faibussowitsch PetscCall(PetscSubcommSetType(red->psubcomm, PETSC_SUBCOMM_CONTIGUOUS));
391e5acf8a4SHong Zhang
3929566063dSJacob Faibussowitsch PetscCall(PetscSubcommSetOptionsPrefix(red->psubcomm, prefix));
3939566063dSJacob Faibussowitsch PetscCall(PetscSubcommSetFromOptions(red->psubcomm));
39475024027SHong Zhang
39575024027SHong Zhang /* create a new PC that processors in each subcomm have copy of */
39675024027SHong Zhang subcomm = PetscSubcommChild(red->psubcomm);
39775024027SHong Zhang
3989566063dSJacob Faibussowitsch PetscCall(KSPCreate(subcomm, &red->ksp));
3993821be0aSBarry Smith PetscCall(KSPSetNestLevel(red->ksp, pc->kspnestlevel));
4009566063dSJacob Faibussowitsch PetscCall(KSPSetErrorIfNotConverged(red->ksp, pc->erroriffailure));
4019566063dSJacob Faibussowitsch PetscCall(PetscObjectIncrementTabLevel((PetscObject)red->ksp, (PetscObject)pc, 1));
4029566063dSJacob Faibussowitsch PetscCall(KSPSetType(red->ksp, KSPPREONLY));
4039566063dSJacob Faibussowitsch PetscCall(KSPGetPC(red->ksp, &red->pc));
4049566063dSJacob Faibussowitsch PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATSEQSBAIJ, &issbaij));
40548a46eb9SPierre Jolivet if (!issbaij) PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATMPISBAIJ, &issbaij));
40608cecb0aSPierre Jolivet if (!issbaij) {
4079566063dSJacob Faibussowitsch PetscCall(PCSetType(red->pc, PCLU));
40808cecb0aSPierre Jolivet } else {
4099566063dSJacob Faibussowitsch PetscCall(PCSetType(red->pc, PCCHOLESKY));
41008cecb0aSPierre Jolivet }
411753b7fb9SBarry Smith if (red->shifttypeset) {
4129566063dSJacob Faibussowitsch PetscCall(PCFactorSetShiftType(red->pc, red->shifttype));
413753b7fb9SBarry Smith red->shifttypeset = PETSC_FALSE;
414753b7fb9SBarry Smith }
4159566063dSJacob Faibussowitsch PetscCall(KSPSetOptionsPrefix(red->ksp, prefix));
4169566063dSJacob Faibussowitsch PetscCall(KSPAppendOptionsPrefix(red->ksp, "redundant_"));
41775024027SHong Zhang }
41883ab6a24SBarry Smith *innerksp = red->ksp;
4193ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
4204b9ad928SBarry Smith }
4214b9ad928SBarry Smith
4224b9ad928SBarry Smith /*@
423f1580f4eSBarry Smith PCRedundantGetKSP - Gets the less parallel `KSP` created by the redundant `PC`.
4244b9ad928SBarry Smith
4254b9ad928SBarry Smith Not Collective
4264b9ad928SBarry Smith
4274b9ad928SBarry Smith Input Parameter:
4284b9ad928SBarry Smith . pc - the preconditioner context
4294b9ad928SBarry Smith
4304b9ad928SBarry Smith Output Parameter:
431f1580f4eSBarry Smith . innerksp - the `KSP` on the smaller set of processes
4324b9ad928SBarry Smith
4334b9ad928SBarry Smith Level: advanced
4344b9ad928SBarry Smith
435562efe2eSBarry Smith .seealso: [](ch_ksp), `PCREDUNDANT`
4364b9ad928SBarry Smith @*/
PCRedundantGetKSP(PC pc,KSP * innerksp)437d71ae5a4SJacob Faibussowitsch PetscErrorCode PCRedundantGetKSP(PC pc, KSP *innerksp)
438d71ae5a4SJacob Faibussowitsch {
4394b9ad928SBarry Smith PetscFunctionBegin;
4400700a824SBarry Smith PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
4414f572ea9SToby Isaac PetscAssertPointer(innerksp, 2);
442cac4c232SBarry Smith PetscUseMethod(pc, "PCRedundantGetKSP_C", (PC, KSP *), (pc, innerksp));
4433ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
4444b9ad928SBarry Smith }
4454b9ad928SBarry Smith
PCRedundantGetOperators_Redundant(PC pc,Mat * mat,Mat * pmat)446d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCRedundantGetOperators_Redundant(PC pc, Mat *mat, Mat *pmat)
447d71ae5a4SJacob Faibussowitsch {
4484b9ad928SBarry Smith PC_Redundant *red = (PC_Redundant *)pc->data;
4494b9ad928SBarry Smith
4504b9ad928SBarry Smith PetscFunctionBegin;
451b3804887SHong Zhang if (mat) *mat = red->pmats;
452b3804887SHong Zhang if (pmat) *pmat = red->pmats;
4533ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
4544b9ad928SBarry Smith }
4554b9ad928SBarry Smith
4564b9ad928SBarry Smith /*@
4577addb90fSBarry Smith PCRedundantGetOperators - gets the sequential linear system matrix and matrix used to construct the preconditioner
4584b9ad928SBarry Smith
4594b9ad928SBarry Smith Not Collective
4604b9ad928SBarry Smith
4614b9ad928SBarry Smith Input Parameter:
4624b9ad928SBarry Smith . pc - the preconditioner context
4634b9ad928SBarry Smith
4644b9ad928SBarry Smith Output Parameters:
4654b9ad928SBarry Smith + mat - the matrix
4667addb90fSBarry Smith - pmat - the (possibly different) matrix used to construct the preconditioner
4674b9ad928SBarry Smith
4684b9ad928SBarry Smith Level: advanced
4694b9ad928SBarry Smith
470562efe2eSBarry Smith .seealso: [](ch_ksp), `PCREDUNDANT`
4714b9ad928SBarry Smith @*/
PCRedundantGetOperators(PC pc,Mat * mat,Mat * pmat)472d71ae5a4SJacob Faibussowitsch PetscErrorCode PCRedundantGetOperators(PC pc, Mat *mat, Mat *pmat)
473d71ae5a4SJacob Faibussowitsch {
4744b9ad928SBarry Smith PetscFunctionBegin;
4750700a824SBarry Smith PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
4764f572ea9SToby Isaac if (mat) PetscAssertPointer(mat, 2);
4774f572ea9SToby Isaac if (pmat) PetscAssertPointer(pmat, 3);
478cac4c232SBarry Smith PetscUseMethod(pc, "PCRedundantGetOperators_C", (PC, Mat *, Mat *), (pc, mat, pmat));
4793ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
4804b9ad928SBarry Smith }
4814b9ad928SBarry Smith
48237a17b4dSBarry Smith /*MC
483f1580f4eSBarry Smith PCREDUNDANT - Runs a `KSP` solver with preconditioner for the entire problem on subgroups of processors
48437a17b4dSBarry Smith
485f1580f4eSBarry Smith Options Database Key:
4869b21d695SBarry Smith . -pc_redundant_number <n> - number of redundant solves, for example if you are using 64 MPI processes and
4879b21d695SBarry Smith use an n of 4 there will be 4 parallel solves each on 16 = 64/4 processes.
48809391456SBarry Smith
48937a17b4dSBarry Smith Level: intermediate
49037a17b4dSBarry Smith
49195452b02SPatrick Sanan Notes:
492f1580f4eSBarry Smith Options for the redundant preconditioners can be set using the options database prefix -redundant_
49383ab6a24SBarry Smith
494f1580f4eSBarry Smith The default `KSP` is preonly and the default `PC` is `PCLU` or `PCCHOLESKY` if Pmat is of type `MATSBAIJ`.
495753b7fb9SBarry Smith
496f1580f4eSBarry Smith `PCFactorSetShiftType()` applied to this `PC` will convey they shift type into the inner `PC` if it is factorization based.
497f1580f4eSBarry Smith
498f1580f4eSBarry Smith Developer Note:
499f1580f4eSBarry Smith `PCSetInitialGuessNonzero()` is not used by this class but likely should be.
5009cfaa89bSBarry Smith
501562efe2eSBarry Smith .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PCRedundantSetScatter()`,
502f1580f4eSBarry Smith `PCRedundantGetKSP()`, `PCRedundantGetOperators()`, `PCRedundantSetNumber()`, `PCREDISTRIBUTE`
50337a17b4dSBarry Smith M*/
50437a17b4dSBarry Smith
PCCreate_Redundant(PC pc)505d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode PCCreate_Redundant(PC pc)
506d71ae5a4SJacob Faibussowitsch {
5074b9ad928SBarry Smith PC_Redundant *red;
50869db28dcSHong Zhang PetscMPIInt size;
5093f457be1SHong Zhang
5104b9ad928SBarry Smith PetscFunctionBegin;
5114dfa11a4SJacob Faibussowitsch PetscCall(PetscNew(&red));
5129566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
5132fa5cd67SKarl Rupp
51469db28dcSHong Zhang red->nsubcomm = size;
5154b9ad928SBarry Smith red->useparallelmat = PETSC_TRUE;
5161fbd8f88SHong Zhang pc->data = (void *)red;
5174b9ad928SBarry Smith
5184b9ad928SBarry Smith pc->ops->apply = PCApply_Redundant;
519d88bfacbSStefano Zampini pc->ops->applytranspose = PCApplyTranspose_Redundant;
5204b9ad928SBarry Smith pc->ops->setup = PCSetUp_Redundant;
5214b9ad928SBarry Smith pc->ops->destroy = PCDestroy_Redundant;
5221ea5a559SBarry Smith pc->ops->reset = PCReset_Redundant;
5234b9ad928SBarry Smith pc->ops->setfromoptions = PCSetFromOptions_Redundant;
5244b9ad928SBarry Smith pc->ops->view = PCView_Redundant;
5252fa5cd67SKarl Rupp
5269566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantSetScatter_C", PCRedundantSetScatter_Redundant));
5279566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantSetNumber_C", PCRedundantSetNumber_Redundant));
5289566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantGetKSP_C", PCRedundantGetKSP_Redundant));
5299566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantGetOperators_C", PCRedundantGetOperators_Redundant));
5309566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFactorSetShiftType_C", PCFactorSetShiftType_Redundant));
5313ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
5324b9ad928SBarry Smith }
533