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