14b9ad928SBarry Smith /* 23f457be1SHong Zhang This file defines a "solve the problem redundantly on each subgroup of processor" preconditioner. 34b9ad928SBarry Smith */ 4af0996ceSBarry Smith #include <petsc/private/pcimpl.h> 5c6db04a5SJed Brown #include <petscksp.h> /*I "petscksp.h" I*/ 64b9ad928SBarry Smith 74b9ad928SBarry Smith typedef struct { 83e065800SHong Zhang KSP ksp; 94b9ad928SBarry Smith PC pc; /* actual preconditioner used on each processor */ 10ce94432eSBarry Smith Vec xsub, ysub; /* vectors of a subcommunicator to hold parallel vectors of PetscObjectComm((PetscObject)pc) */ 113f457be1SHong Zhang Vec xdup, ydup; /* parallel vector that congregates xsub or ysub facilitating vector scattering */ 12b3804887SHong Zhang Mat pmats; /* matrix and optional preconditioner matrix belong to a subcommunicator */ 133f457be1SHong Zhang VecScatter scatterin, scatterout; /* scatter used to move all values to each processor group (subcommunicator) */ 14ace3abfcSBarry Smith PetscBool useparallelmat; 15c540e29cSHong Zhang PetscSubcomm psubcomm; 161fbd8f88SHong Zhang PetscInt nsubcomm; /* num of data structure PetscSubcomm */ 17753b7fb9SBarry Smith PetscBool shifttypeset; 18753b7fb9SBarry Smith MatFactorShiftType shifttype; 194b9ad928SBarry Smith } PC_Redundant; 204b9ad928SBarry Smith 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 37d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCView_Redundant(PC pc, PetscViewer viewer) 38d71ae5a4SJacob Faibussowitsch { 394b9ad928SBarry Smith PC_Redundant *red = (PC_Redundant *)pc->data; 40ace3abfcSBarry Smith PetscBool iascii, isstring; 4103ccd0b4SBarry Smith PetscViewer subviewer; 424b9ad928SBarry Smith 434b9ad928SBarry Smith PetscFunctionBegin; 449566063dSJacob Faibussowitsch PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 459566063dSJacob Faibussowitsch PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring)); 4632077d6dSBarry Smith if (iascii) { 4703ccd0b4SBarry Smith if (!red->psubcomm) { 489566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " Not yet setup\n")); 4903ccd0b4SBarry Smith } else { 5063a3b9bcSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " First (color=0) of %" PetscInt_FMT " PCs follows\n", red->nsubcomm)); 519566063dSJacob Faibussowitsch PetscCall(PetscViewerGetSubViewer(viewer, ((PetscObject)red->pc)->comm, &subviewer)); 52f5dd71faSBarry Smith if (!red->psubcomm->color) { /* only view first redundant pc */ 539566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPushTab(subviewer)); 549566063dSJacob Faibussowitsch PetscCall(KSPView(red->ksp, subviewer)); 559566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPopTab(subviewer)); 564b9ad928SBarry Smith } 579566063dSJacob Faibussowitsch PetscCall(PetscViewerRestoreSubViewer(viewer, ((PetscObject)red->pc)->comm, &subviewer)); 584b9ad928SBarry Smith } 5903ccd0b4SBarry Smith } else if (isstring) { 609566063dSJacob Faibussowitsch PetscCall(PetscViewerStringSPrintf(viewer, " Redundant solver preconditioner")); 614b9ad928SBarry Smith } 623ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 634b9ad928SBarry Smith } 644b9ad928SBarry Smith 6519b3b6edSHong Zhang #include <../src/mat/impls/aij/mpi/mpiaij.h> 66d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCSetUp_Redundant(PC pc) 67d71ae5a4SJacob Faibussowitsch { 684b9ad928SBarry Smith PC_Redundant *red = (PC_Redundant *)pc->data; 691b81debcSHong Zhang PetscInt mstart, mend, mlocal, M; 7013f74950SBarry Smith PetscMPIInt size; 71ce94432eSBarry Smith MPI_Comm comm, subcomm; 72ddc54837SHong Zhang Vec x; 733f457be1SHong Zhang 744b9ad928SBarry Smith PetscFunctionBegin; 759566063dSJacob Faibussowitsch PetscCall(PetscObjectGetComm((PetscObject)pc, &comm)); 76ddc54837SHong Zhang 77ddc54837SHong Zhang /* if pmatrix set by user is sequential then we do not need to gather the parallel matrix */ 789566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_size(comm, &size)); 79ddc54837SHong Zhang if (size == 1) red->useparallelmat = PETSC_FALSE; 801fbd8f88SHong Zhang 814b9ad928SBarry Smith if (!pc->setupcalled) { 821b81debcSHong Zhang PetscInt mloc_sub; 8375024027SHong Zhang if (!red->psubcomm) { /* create red->psubcomm, new ksp and pc over subcomm */ 8475024027SHong Zhang KSP ksp; 859566063dSJacob Faibussowitsch PetscCall(PCRedundantGetKSP(pc, &ksp)); 861b81debcSHong Zhang } 8775024027SHong Zhang subcomm = PetscSubcommChild(red->psubcomm); 881fbd8f88SHong Zhang 891b81debcSHong Zhang if (red->useparallelmat) { 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. 1168b96b0d2SHong Zhang xdup concatenates all xsub's contigously to form a mpi vector over dupcomm (see PetscSubcommCreate_interlaced()) 1178b96b0d2SHong Zhang ydup concatenates all ysub and has empty local arrays because ysub's arrays will be place into it. 118ce94432eSBarry Smith Note: we use communicator dupcomm, not PetscObjectComm((PetscObject)pc)! */ 1199566063dSJacob Faibussowitsch PetscCall(MatGetLocalSize(red->pmats, &mloc_sub, NULL)); 1209566063dSJacob Faibussowitsch PetscCall(VecCreateMPI(PetscSubcommContiguousParent(red->psubcomm), mloc_sub, PETSC_DECIDE, &red->xdup)); 1219566063dSJacob Faibussowitsch PetscCall(VecCreateMPIWithArray(PetscSubcommContiguousParent(red->psubcomm), 1, mloc_sub, PETSC_DECIDE, NULL, &red->ydup)); 1223f457be1SHong Zhang 123f68be91cSHong Zhang /* create vecscatters */ 124f68be91cSHong Zhang if (!red->scatterin) { /* efficiency of scatterin is independent from psubcomm_type! */ 1253f457be1SHong Zhang IS is1, is2; 1263f457be1SHong Zhang PetscInt *idx1, *idx2, i, j, k; 12745fc02eaSBarry Smith 1289566063dSJacob Faibussowitsch PetscCall(MatCreateVecs(pc->pmat, &x, NULL)); 1299566063dSJacob Faibussowitsch PetscCall(VecGetSize(x, &M)); 1309566063dSJacob Faibussowitsch PetscCall(VecGetOwnershipRange(x, &mstart, &mend)); 1311b81debcSHong Zhang mlocal = mend - mstart; 1329566063dSJacob Faibussowitsch PetscCall(PetscMalloc2(red->psubcomm->n * mlocal, &idx1, red->psubcomm->n * mlocal, &idx2)); 1333f457be1SHong Zhang j = 0; 1341fbd8f88SHong Zhang for (k = 0; k < red->psubcomm->n; k++) { 1353f457be1SHong Zhang for (i = mstart; i < mend; i++) { 1363f457be1SHong Zhang idx1[j] = i; 137ddc54837SHong Zhang idx2[j++] = i + M * k; 1383f457be1SHong Zhang } 1393f457be1SHong Zhang } 1409566063dSJacob Faibussowitsch PetscCall(ISCreateGeneral(comm, red->psubcomm->n * mlocal, idx1, PETSC_COPY_VALUES, &is1)); 1419566063dSJacob Faibussowitsch PetscCall(ISCreateGeneral(comm, red->psubcomm->n * mlocal, idx2, PETSC_COPY_VALUES, &is2)); 1429566063dSJacob Faibussowitsch PetscCall(VecScatterCreate(x, is1, red->xdup, is2, &red->scatterin)); 1439566063dSJacob Faibussowitsch PetscCall(ISDestroy(&is1)); 1449566063dSJacob Faibussowitsch PetscCall(ISDestroy(&is2)); 1453f457be1SHong Zhang 1466909748bSHong Zhang /* Impl below is good for PETSC_SUBCOMM_INTERLACED (no inter-process communication) and PETSC_SUBCOMM_CONTIGUOUS (communication within subcomm) */ 1479566063dSJacob Faibussowitsch PetscCall(ISCreateStride(comm, mlocal, mstart + red->psubcomm->color * M, 1, &is1)); 1489566063dSJacob Faibussowitsch PetscCall(ISCreateStride(comm, mlocal, mstart, 1, &is2)); 1499566063dSJacob Faibussowitsch PetscCall(VecScatterCreate(red->xdup, is1, x, is2, &red->scatterout)); 1509566063dSJacob Faibussowitsch PetscCall(ISDestroy(&is1)); 1519566063dSJacob Faibussowitsch PetscCall(ISDestroy(&is2)); 1529566063dSJacob Faibussowitsch PetscCall(PetscFree2(idx1, idx2)); 1539566063dSJacob Faibussowitsch PetscCall(VecDestroy(&x)); 1541b81debcSHong Zhang } 155ab661555SHong Zhang } else { /* !red->useparallelmat */ 1569566063dSJacob Faibussowitsch PetscCall(KSPSetOperators(red->ksp, pc->mat, pc->pmat)); 1571b81debcSHong Zhang } 158ab661555SHong Zhang } else { /* pc->setupcalled */ 1594b9ad928SBarry Smith if (red->useparallelmat) { 160ab661555SHong Zhang MatReuse reuse; 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)); 1779566063dSJacob Faibussowitsch PetscCall(KSPSetUp(red->ksp)); 1789767e1d8SStefano Zampini 1799767e1d8SStefano Zampini /* Detect failure */ 1809767e1d8SStefano Zampini KSPConvergedReason redreason; 1819767e1d8SStefano Zampini PetscCall(KSPGetConvergedReason(red->ksp, &redreason)); 1829767e1d8SStefano Zampini if (redreason) pc->failedreason = PC_SUBPC_ERROR; 1833ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1844b9ad928SBarry Smith } 1854b9ad928SBarry Smith 186d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCApply_Redundant(PC pc, Vec x, Vec y) 187d71ae5a4SJacob Faibussowitsch { 1884b9ad928SBarry Smith PC_Redundant *red = (PC_Redundant *)pc->data; 1893f457be1SHong Zhang PetscScalar *array; 1904b9ad928SBarry Smith 1914b9ad928SBarry Smith PetscFunctionBegin; 192ddc54837SHong Zhang if (!red->useparallelmat) { 1939566063dSJacob Faibussowitsch PetscCall(KSPSolve(red->ksp, x, y)); 1949566063dSJacob Faibussowitsch PetscCall(KSPCheckSolve(red->ksp, pc, y)); 1953ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 196ddc54837SHong Zhang } 197ddc54837SHong Zhang 1983f457be1SHong Zhang /* scatter x to xdup */ 1999566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(red->scatterin, x, red->xdup, INSERT_VALUES, SCATTER_FORWARD)); 2009566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(red->scatterin, x, red->xdup, INSERT_VALUES, SCATTER_FORWARD)); 2013f457be1SHong Zhang 2023f457be1SHong Zhang /* place xdup's local array into xsub */ 2039566063dSJacob Faibussowitsch PetscCall(VecGetArray(red->xdup, &array)); 2049566063dSJacob Faibussowitsch PetscCall(VecPlaceArray(red->xsub, (const PetscScalar *)array)); 2054b9ad928SBarry Smith 2064b9ad928SBarry Smith /* apply preconditioner on each processor */ 2079566063dSJacob Faibussowitsch PetscCall(KSPSolve(red->ksp, red->xsub, red->ysub)); 2089566063dSJacob Faibussowitsch PetscCall(KSPCheckSolve(red->ksp, pc, red->ysub)); 2099566063dSJacob Faibussowitsch PetscCall(VecResetArray(red->xsub)); 2109566063dSJacob Faibussowitsch PetscCall(VecRestoreArray(red->xdup, &array)); 2114b9ad928SBarry Smith 2123f457be1SHong Zhang /* place ysub's local array into ydup */ 2139566063dSJacob Faibussowitsch PetscCall(VecGetArray(red->ysub, &array)); 2149566063dSJacob Faibussowitsch PetscCall(VecPlaceArray(red->ydup, (const PetscScalar *)array)); 2153f457be1SHong Zhang 2163f457be1SHong Zhang /* scatter ydup to y */ 2179566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(red->scatterout, red->ydup, y, INSERT_VALUES, SCATTER_FORWARD)); 2189566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(red->scatterout, red->ydup, y, INSERT_VALUES, SCATTER_FORWARD)); 2199566063dSJacob Faibussowitsch PetscCall(VecResetArray(red->ydup)); 2209566063dSJacob Faibussowitsch PetscCall(VecRestoreArray(red->ysub, &array)); 2213ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 2224b9ad928SBarry Smith } 2234b9ad928SBarry Smith 224d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCApplyTranspose_Redundant(PC pc, Vec x, Vec y) 225d71ae5a4SJacob Faibussowitsch { 226d88bfacbSStefano Zampini PC_Redundant *red = (PC_Redundant *)pc->data; 227d88bfacbSStefano Zampini PetscScalar *array; 228d88bfacbSStefano Zampini 229d88bfacbSStefano Zampini PetscFunctionBegin; 230d88bfacbSStefano Zampini if (!red->useparallelmat) { 2319566063dSJacob Faibussowitsch PetscCall(KSPSolveTranspose(red->ksp, x, y)); 2329566063dSJacob Faibussowitsch PetscCall(KSPCheckSolve(red->ksp, pc, y)); 2333ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 234d88bfacbSStefano Zampini } 235d88bfacbSStefano Zampini 236d88bfacbSStefano Zampini /* scatter x to xdup */ 2379566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(red->scatterin, x, red->xdup, INSERT_VALUES, SCATTER_FORWARD)); 2389566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(red->scatterin, x, red->xdup, INSERT_VALUES, SCATTER_FORWARD)); 239d88bfacbSStefano Zampini 240d88bfacbSStefano Zampini /* place xdup's local array into xsub */ 2419566063dSJacob Faibussowitsch PetscCall(VecGetArray(red->xdup, &array)); 2429566063dSJacob Faibussowitsch PetscCall(VecPlaceArray(red->xsub, (const PetscScalar *)array)); 243d88bfacbSStefano Zampini 244d88bfacbSStefano Zampini /* apply preconditioner on each processor */ 2459566063dSJacob Faibussowitsch PetscCall(KSPSolveTranspose(red->ksp, red->xsub, red->ysub)); 2469566063dSJacob Faibussowitsch PetscCall(KSPCheckSolve(red->ksp, pc, red->ysub)); 2479566063dSJacob Faibussowitsch PetscCall(VecResetArray(red->xsub)); 2489566063dSJacob Faibussowitsch PetscCall(VecRestoreArray(red->xdup, &array)); 249d88bfacbSStefano Zampini 250d88bfacbSStefano Zampini /* place ysub's local array into ydup */ 2519566063dSJacob Faibussowitsch PetscCall(VecGetArray(red->ysub, &array)); 2529566063dSJacob Faibussowitsch PetscCall(VecPlaceArray(red->ydup, (const PetscScalar *)array)); 253d88bfacbSStefano Zampini 254d88bfacbSStefano Zampini /* scatter ydup to y */ 2559566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(red->scatterout, red->ydup, y, INSERT_VALUES, SCATTER_FORWARD)); 2569566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(red->scatterout, red->ydup, y, INSERT_VALUES, SCATTER_FORWARD)); 2579566063dSJacob Faibussowitsch PetscCall(VecResetArray(red->ydup)); 2589566063dSJacob Faibussowitsch PetscCall(VecRestoreArray(red->ysub, &array)); 2593ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 260d88bfacbSStefano Zampini } 261d88bfacbSStefano Zampini 262d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCReset_Redundant(PC pc) 263d71ae5a4SJacob Faibussowitsch { 2644b9ad928SBarry Smith PC_Redundant *red = (PC_Redundant *)pc->data; 2654b9ad928SBarry Smith 2664b9ad928SBarry Smith PetscFunctionBegin; 2671b81debcSHong Zhang if (red->useparallelmat) { 2689566063dSJacob Faibussowitsch PetscCall(VecScatterDestroy(&red->scatterin)); 2699566063dSJacob Faibussowitsch PetscCall(VecScatterDestroy(&red->scatterout)); 2709566063dSJacob Faibussowitsch PetscCall(VecDestroy(&red->ysub)); 2719566063dSJacob Faibussowitsch PetscCall(VecDestroy(&red->xsub)); 2729566063dSJacob Faibussowitsch PetscCall(VecDestroy(&red->xdup)); 2739566063dSJacob Faibussowitsch PetscCall(VecDestroy(&red->ydup)); 2741b81debcSHong Zhang } 2759566063dSJacob Faibussowitsch PetscCall(MatDestroy(&red->pmats)); 2769566063dSJacob Faibussowitsch PetscCall(KSPReset(red->ksp)); 2773ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 2781ea5a559SBarry Smith } 2791ea5a559SBarry Smith 280d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCDestroy_Redundant(PC pc) 281d71ae5a4SJacob Faibussowitsch { 2821ea5a559SBarry Smith PC_Redundant *red = (PC_Redundant *)pc->data; 2831ea5a559SBarry Smith 2841ea5a559SBarry Smith PetscFunctionBegin; 2859566063dSJacob Faibussowitsch PetscCall(PCReset_Redundant(pc)); 2869566063dSJacob Faibussowitsch PetscCall(KSPDestroy(&red->ksp)); 2879566063dSJacob Faibussowitsch PetscCall(PetscSubcommDestroy(&red->psubcomm)); 2882e956fe4SStefano Zampini PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantSetScatter_C", NULL)); 2892e956fe4SStefano Zampini PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantSetNumber_C", NULL)); 2902e956fe4SStefano Zampini PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantGetKSP_C", NULL)); 2912e956fe4SStefano Zampini PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantGetOperators_C", NULL)); 2922e956fe4SStefano Zampini PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFactorSetShiftType_C", NULL)); 2939566063dSJacob Faibussowitsch PetscCall(PetscFree(pc->data)); 2943ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 2954b9ad928SBarry Smith } 2964b9ad928SBarry Smith 297d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCSetFromOptions_Redundant(PC pc, PetscOptionItems *PetscOptionsObject) 298d71ae5a4SJacob Faibussowitsch { 299a98ce0f4SHong Zhang PC_Redundant *red = (PC_Redundant *)pc->data; 300a98ce0f4SHong Zhang 3014b9ad928SBarry Smith PetscFunctionBegin; 302d0609cedSBarry Smith PetscOptionsHeadBegin(PetscOptionsObject, "Redundant options"); 3039566063dSJacob Faibussowitsch PetscCall(PetscOptionsInt("-pc_redundant_number", "Number of redundant pc", "PCRedundantSetNumber", red->nsubcomm, &red->nsubcomm, NULL)); 304d0609cedSBarry Smith PetscOptionsHeadEnd(); 3053ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 3064b9ad928SBarry Smith } 3074b9ad928SBarry Smith 308d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCRedundantSetNumber_Redundant(PC pc, PetscInt nreds) 309d71ae5a4SJacob Faibussowitsch { 31009a6bc64SHong Zhang PC_Redundant *red = (PC_Redundant *)pc->data; 31109a6bc64SHong Zhang 31209a6bc64SHong Zhang PetscFunctionBegin; 31309a6bc64SHong Zhang red->nsubcomm = nreds; 3143ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 31509a6bc64SHong Zhang } 31609a6bc64SHong Zhang 31709a6bc64SHong Zhang /*@ 31809a6bc64SHong Zhang PCRedundantSetNumber - Sets the number of redundant preconditioner contexts. 31909a6bc64SHong Zhang 320c3339decSBarry Smith Logically Collective 32109a6bc64SHong Zhang 32209a6bc64SHong Zhang Input Parameters: 32309a6bc64SHong Zhang + pc - the preconditioner context 3249b21d695SBarry Smith - nredundant - number of redundant preconditioner contexts; for example if you are using 64 MPI processes and 3259b21d695SBarry Smith use an nredundant of 4 there will be 4 parallel solves each on 16 = 64/4 processes. 32609a6bc64SHong Zhang 32709a6bc64SHong Zhang Level: advanced 32809a6bc64SHong Zhang 329*562efe2eSBarry Smith .seealso: [](ch_ksp), `PCREDUNDANT` 33009a6bc64SHong Zhang @*/ 331d71ae5a4SJacob Faibussowitsch PetscErrorCode PCRedundantSetNumber(PC pc, PetscInt nredundant) 332d71ae5a4SJacob Faibussowitsch { 33309a6bc64SHong Zhang PetscFunctionBegin; 3340700a824SBarry Smith PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 33563a3b9bcSJacob Faibussowitsch PetscCheck(nredundant > 0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "num of redundant pc %" PetscInt_FMT " must be positive", nredundant); 336cac4c232SBarry Smith PetscTryMethod(pc, "PCRedundantSetNumber_C", (PC, PetscInt), (pc, nredundant)); 3373ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 33809a6bc64SHong Zhang } 33909a6bc64SHong Zhang 340d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCRedundantSetScatter_Redundant(PC pc, VecScatter in, VecScatter out) 341d71ae5a4SJacob Faibussowitsch { 3424b9ad928SBarry Smith PC_Redundant *red = (PC_Redundant *)pc->data; 3434b9ad928SBarry Smith 3444b9ad928SBarry Smith PetscFunctionBegin; 3459566063dSJacob Faibussowitsch PetscCall(PetscObjectReference((PetscObject)in)); 3469566063dSJacob Faibussowitsch PetscCall(VecScatterDestroy(&red->scatterin)); 3472fa5cd67SKarl Rupp 348c3122656SLisandro Dalcin red->scatterin = in; 3492fa5cd67SKarl Rupp 3509566063dSJacob Faibussowitsch PetscCall(PetscObjectReference((PetscObject)out)); 3519566063dSJacob Faibussowitsch PetscCall(VecScatterDestroy(&red->scatterout)); 352c3122656SLisandro Dalcin red->scatterout = out; 3533ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 3544b9ad928SBarry Smith } 3554b9ad928SBarry Smith 3564b9ad928SBarry Smith /*@ 3574b9ad928SBarry Smith PCRedundantSetScatter - Sets the scatter used to copy values into the 3584b9ad928SBarry Smith redundant local solve and the scatter to move them back into the global 3594b9ad928SBarry Smith vector. 3604b9ad928SBarry Smith 361c3339decSBarry Smith Logically Collective 3624b9ad928SBarry Smith 3634b9ad928SBarry Smith Input Parameters: 3644b9ad928SBarry Smith + pc - the preconditioner context 3654b9ad928SBarry Smith . in - the scatter to move the values in 3664b9ad928SBarry Smith - out - the scatter to move them out 3674b9ad928SBarry Smith 3684b9ad928SBarry Smith Level: advanced 3694b9ad928SBarry Smith 370*562efe2eSBarry Smith .seealso: [](ch_ksp), `PCREDUNDANT` 3714b9ad928SBarry Smith @*/ 372d71ae5a4SJacob Faibussowitsch PetscErrorCode PCRedundantSetScatter(PC pc, VecScatter in, VecScatter out) 373d71ae5a4SJacob Faibussowitsch { 3744b9ad928SBarry Smith PetscFunctionBegin; 3750700a824SBarry Smith PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 37697929ea7SJunchao Zhang PetscValidHeaderSpecific(in, PETSCSF_CLASSID, 2); 37797929ea7SJunchao Zhang PetscValidHeaderSpecific(out, PETSCSF_CLASSID, 3); 378cac4c232SBarry Smith PetscTryMethod(pc, "PCRedundantSetScatter_C", (PC, VecScatter, VecScatter), (pc, in, out)); 3793ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 3804b9ad928SBarry Smith } 3814b9ad928SBarry Smith 382d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCRedundantGetKSP_Redundant(PC pc, KSP *innerksp) 383d71ae5a4SJacob Faibussowitsch { 3844b9ad928SBarry Smith PC_Redundant *red = (PC_Redundant *)pc->data; 38575024027SHong Zhang MPI_Comm comm, subcomm; 38675024027SHong Zhang const char *prefix; 38708cecb0aSPierre Jolivet PetscBool issbaij; 3884b9ad928SBarry Smith 3894b9ad928SBarry Smith PetscFunctionBegin; 39075024027SHong Zhang if (!red->psubcomm) { 3919566063dSJacob Faibussowitsch PetscCall(PCGetOptionsPrefix(pc, &prefix)); 392e5acf8a4SHong Zhang 3939566063dSJacob Faibussowitsch PetscCall(PetscObjectGetComm((PetscObject)pc, &comm)); 3949566063dSJacob Faibussowitsch PetscCall(PetscSubcommCreate(comm, &red->psubcomm)); 3959566063dSJacob Faibussowitsch PetscCall(PetscSubcommSetNumber(red->psubcomm, red->nsubcomm)); 3969566063dSJacob Faibussowitsch PetscCall(PetscSubcommSetType(red->psubcomm, PETSC_SUBCOMM_CONTIGUOUS)); 397e5acf8a4SHong Zhang 3989566063dSJacob Faibussowitsch PetscCall(PetscSubcommSetOptionsPrefix(red->psubcomm, prefix)); 3999566063dSJacob Faibussowitsch PetscCall(PetscSubcommSetFromOptions(red->psubcomm)); 40075024027SHong Zhang 40175024027SHong Zhang /* create a new PC that processors in each subcomm have copy of */ 40275024027SHong Zhang subcomm = PetscSubcommChild(red->psubcomm); 40375024027SHong Zhang 4049566063dSJacob Faibussowitsch PetscCall(KSPCreate(subcomm, &red->ksp)); 4053821be0aSBarry Smith PetscCall(KSPSetNestLevel(red->ksp, pc->kspnestlevel)); 4069566063dSJacob Faibussowitsch PetscCall(KSPSetErrorIfNotConverged(red->ksp, pc->erroriffailure)); 4079566063dSJacob Faibussowitsch PetscCall(PetscObjectIncrementTabLevel((PetscObject)red->ksp, (PetscObject)pc, 1)); 4089566063dSJacob Faibussowitsch PetscCall(KSPSetType(red->ksp, KSPPREONLY)); 4099566063dSJacob Faibussowitsch PetscCall(KSPGetPC(red->ksp, &red->pc)); 4109566063dSJacob Faibussowitsch PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATSEQSBAIJ, &issbaij)); 41148a46eb9SPierre Jolivet if (!issbaij) PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATMPISBAIJ, &issbaij)); 41208cecb0aSPierre Jolivet if (!issbaij) { 4139566063dSJacob Faibussowitsch PetscCall(PCSetType(red->pc, PCLU)); 41408cecb0aSPierre Jolivet } else { 4159566063dSJacob Faibussowitsch PetscCall(PCSetType(red->pc, PCCHOLESKY)); 41608cecb0aSPierre Jolivet } 417753b7fb9SBarry Smith if (red->shifttypeset) { 4189566063dSJacob Faibussowitsch PetscCall(PCFactorSetShiftType(red->pc, red->shifttype)); 419753b7fb9SBarry Smith red->shifttypeset = PETSC_FALSE; 420753b7fb9SBarry Smith } 4219566063dSJacob Faibussowitsch PetscCall(KSPSetOptionsPrefix(red->ksp, prefix)); 4229566063dSJacob Faibussowitsch PetscCall(KSPAppendOptionsPrefix(red->ksp, "redundant_")); 42375024027SHong Zhang } 42483ab6a24SBarry Smith *innerksp = red->ksp; 4253ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 4264b9ad928SBarry Smith } 4274b9ad928SBarry Smith 4284b9ad928SBarry Smith /*@ 429f1580f4eSBarry Smith PCRedundantGetKSP - Gets the less parallel `KSP` created by the redundant `PC`. 4304b9ad928SBarry Smith 4314b9ad928SBarry Smith Not Collective 4324b9ad928SBarry Smith 4334b9ad928SBarry Smith Input Parameter: 4344b9ad928SBarry Smith . pc - the preconditioner context 4354b9ad928SBarry Smith 4364b9ad928SBarry Smith Output Parameter: 437f1580f4eSBarry Smith . innerksp - the `KSP` on the smaller set of processes 4384b9ad928SBarry Smith 4394b9ad928SBarry Smith Level: advanced 4404b9ad928SBarry Smith 441*562efe2eSBarry Smith .seealso: [](ch_ksp), `PCREDUNDANT` 4424b9ad928SBarry Smith @*/ 443d71ae5a4SJacob Faibussowitsch PetscErrorCode PCRedundantGetKSP(PC pc, KSP *innerksp) 444d71ae5a4SJacob Faibussowitsch { 4454b9ad928SBarry Smith PetscFunctionBegin; 4460700a824SBarry Smith PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 4474f572ea9SToby Isaac PetscAssertPointer(innerksp, 2); 448cac4c232SBarry Smith PetscUseMethod(pc, "PCRedundantGetKSP_C", (PC, KSP *), (pc, innerksp)); 4493ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 4504b9ad928SBarry Smith } 4514b9ad928SBarry Smith 452d71ae5a4SJacob Faibussowitsch static PetscErrorCode PCRedundantGetOperators_Redundant(PC pc, Mat *mat, Mat *pmat) 453d71ae5a4SJacob Faibussowitsch { 4544b9ad928SBarry Smith PC_Redundant *red = (PC_Redundant *)pc->data; 4554b9ad928SBarry Smith 4564b9ad928SBarry Smith PetscFunctionBegin; 457b3804887SHong Zhang if (mat) *mat = red->pmats; 458b3804887SHong Zhang if (pmat) *pmat = red->pmats; 4593ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 4604b9ad928SBarry Smith } 4614b9ad928SBarry Smith 4624b9ad928SBarry Smith /*@ 4634b9ad928SBarry Smith PCRedundantGetOperators - gets the sequential matrix and preconditioner matrix 4644b9ad928SBarry Smith 4654b9ad928SBarry Smith Not Collective 4664b9ad928SBarry Smith 4674b9ad928SBarry Smith Input Parameter: 4684b9ad928SBarry Smith . pc - the preconditioner context 4694b9ad928SBarry Smith 4704b9ad928SBarry Smith Output Parameters: 4714b9ad928SBarry Smith + mat - the matrix 4724b9ad928SBarry Smith - pmat - the (possibly different) preconditioner matrix 4734b9ad928SBarry Smith 4744b9ad928SBarry Smith Level: advanced 4754b9ad928SBarry Smith 476*562efe2eSBarry Smith .seealso: [](ch_ksp), `PCREDUNDANT` 4774b9ad928SBarry Smith @*/ 478d71ae5a4SJacob Faibussowitsch PetscErrorCode PCRedundantGetOperators(PC pc, Mat *mat, Mat *pmat) 479d71ae5a4SJacob Faibussowitsch { 4804b9ad928SBarry Smith PetscFunctionBegin; 4810700a824SBarry Smith PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 4824f572ea9SToby Isaac if (mat) PetscAssertPointer(mat, 2); 4834f572ea9SToby Isaac if (pmat) PetscAssertPointer(pmat, 3); 484cac4c232SBarry Smith PetscUseMethod(pc, "PCRedundantGetOperators_C", (PC, Mat *, Mat *), (pc, mat, pmat)); 4853ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 4864b9ad928SBarry Smith } 4874b9ad928SBarry Smith 48837a17b4dSBarry Smith /*MC 489f1580f4eSBarry Smith PCREDUNDANT - Runs a `KSP` solver with preconditioner for the entire problem on subgroups of processors 49037a17b4dSBarry Smith 491f1580f4eSBarry Smith Options Database Key: 4929b21d695SBarry Smith . -pc_redundant_number <n> - number of redundant solves, for example if you are using 64 MPI processes and 4939b21d695SBarry Smith use an n of 4 there will be 4 parallel solves each on 16 = 64/4 processes. 49409391456SBarry Smith 49537a17b4dSBarry Smith Level: intermediate 49637a17b4dSBarry Smith 49795452b02SPatrick Sanan Notes: 498f1580f4eSBarry Smith Options for the redundant preconditioners can be set using the options database prefix -redundant_ 49983ab6a24SBarry Smith 500f1580f4eSBarry Smith The default `KSP` is preonly and the default `PC` is `PCLU` or `PCCHOLESKY` if Pmat is of type `MATSBAIJ`. 501753b7fb9SBarry Smith 502f1580f4eSBarry Smith `PCFactorSetShiftType()` applied to this `PC` will convey they shift type into the inner `PC` if it is factorization based. 503f1580f4eSBarry Smith 504f1580f4eSBarry Smith Developer Note: 505f1580f4eSBarry Smith `PCSetInitialGuessNonzero()` is not used by this class but likely should be. 5069cfaa89bSBarry Smith 507*562efe2eSBarry Smith .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PCRedundantSetScatter()`, 508f1580f4eSBarry Smith `PCRedundantGetKSP()`, `PCRedundantGetOperators()`, `PCRedundantSetNumber()`, `PCREDISTRIBUTE` 50937a17b4dSBarry Smith M*/ 51037a17b4dSBarry Smith 511d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode PCCreate_Redundant(PC pc) 512d71ae5a4SJacob Faibussowitsch { 5134b9ad928SBarry Smith PC_Redundant *red; 51469db28dcSHong Zhang PetscMPIInt size; 5153f457be1SHong Zhang 5164b9ad928SBarry Smith PetscFunctionBegin; 5174dfa11a4SJacob Faibussowitsch PetscCall(PetscNew(&red)); 5189566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size)); 5192fa5cd67SKarl Rupp 52069db28dcSHong Zhang red->nsubcomm = size; 5214b9ad928SBarry Smith red->useparallelmat = PETSC_TRUE; 5221fbd8f88SHong Zhang pc->data = (void *)red; 5234b9ad928SBarry Smith 5244b9ad928SBarry Smith pc->ops->apply = PCApply_Redundant; 525d88bfacbSStefano Zampini pc->ops->applytranspose = PCApplyTranspose_Redundant; 5264b9ad928SBarry Smith pc->ops->setup = PCSetUp_Redundant; 5274b9ad928SBarry Smith pc->ops->destroy = PCDestroy_Redundant; 5281ea5a559SBarry Smith pc->ops->reset = PCReset_Redundant; 5294b9ad928SBarry Smith pc->ops->setfromoptions = PCSetFromOptions_Redundant; 5304b9ad928SBarry Smith pc->ops->view = PCView_Redundant; 5312fa5cd67SKarl Rupp 5329566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantSetScatter_C", PCRedundantSetScatter_Redundant)); 5339566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantSetNumber_C", PCRedundantSetNumber_Redundant)); 5349566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantGetKSP_C", PCRedundantGetKSP_Redundant)); 5359566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCRedundantGetOperators_C", PCRedundantGetOperators_Redundant)); 5369566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFactorSetShiftType_C", PCFactorSetShiftType_Redundant)); 5373ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 5384b9ad928SBarry Smith } 539