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 22753b7fb9SBarry Smith PetscErrorCode PCFactorSetShiftType_Redundant(PC pc,MatFactorShiftType shifttype) 23753b7fb9SBarry Smith { 24753b7fb9SBarry Smith PC_Redundant *red = (PC_Redundant*)pc->data; 25753b7fb9SBarry Smith 26753b7fb9SBarry Smith PetscFunctionBegin; 27753b7fb9SBarry Smith if (red->ksp) { 28753b7fb9SBarry Smith PC pc; 299566063dSJacob Faibussowitsch PetscCall(KSPGetPC(red->ksp,&pc)); 309566063dSJacob Faibussowitsch PetscCall(PCFactorSetShiftType(pc,shifttype)); 31753b7fb9SBarry Smith } else { 32753b7fb9SBarry Smith red->shifttypeset = PETSC_TRUE; 33753b7fb9SBarry Smith red->shifttype = shifttype; 34753b7fb9SBarry Smith } 35753b7fb9SBarry Smith PetscFunctionReturn(0); 36753b7fb9SBarry Smith } 37753b7fb9SBarry Smith 386849ba73SBarry Smith static PetscErrorCode PCView_Redundant(PC pc,PetscViewer viewer) 394b9ad928SBarry Smith { 404b9ad928SBarry Smith PC_Redundant *red = (PC_Redundant*)pc->data; 41ace3abfcSBarry Smith PetscBool iascii,isstring; 4203ccd0b4SBarry Smith PetscViewer subviewer; 434b9ad928SBarry Smith 444b9ad928SBarry Smith PetscFunctionBegin; 459566063dSJacob Faibussowitsch PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 469566063dSJacob Faibussowitsch PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring)); 4732077d6dSBarry Smith if (iascii) { 4803ccd0b4SBarry Smith if (!red->psubcomm) { 499566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer," Not yet setup\n")); 5003ccd0b4SBarry Smith } else { 5163a3b9bcSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer," First (color=0) of %" PetscInt_FMT " PCs follows\n",red->nsubcomm)); 529566063dSJacob Faibussowitsch PetscCall(PetscViewerGetSubViewer(viewer,((PetscObject)red->pc)->comm,&subviewer)); 53f5dd71faSBarry Smith if (!red->psubcomm->color) { /* only view first redundant pc */ 549566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPushTab(subviewer)); 559566063dSJacob Faibussowitsch PetscCall(KSPView(red->ksp,subviewer)); 569566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPopTab(subviewer)); 574b9ad928SBarry Smith } 589566063dSJacob Faibussowitsch PetscCall(PetscViewerRestoreSubViewer(viewer,((PetscObject)red->pc)->comm,&subviewer)); 594b9ad928SBarry Smith } 6003ccd0b4SBarry Smith } else if (isstring) { 619566063dSJacob Faibussowitsch PetscCall(PetscViewerStringSPrintf(viewer," Redundant solver preconditioner")); 624b9ad928SBarry Smith } 634b9ad928SBarry Smith PetscFunctionReturn(0); 644b9ad928SBarry Smith } 654b9ad928SBarry Smith 6619b3b6edSHong Zhang #include <../src/mat/impls/aij/mpi/mpiaij.h> 676849ba73SBarry Smith static PetscErrorCode PCSetUp_Redundant(PC pc) 684b9ad928SBarry Smith { 694b9ad928SBarry Smith PC_Redundant *red = (PC_Redundant*)pc->data; 701b81debcSHong Zhang PetscInt mstart,mend,mlocal,M; 7113f74950SBarry Smith PetscMPIInt size; 72ce94432eSBarry Smith MPI_Comm comm,subcomm; 73ddc54837SHong Zhang Vec x; 743f457be1SHong Zhang 754b9ad928SBarry Smith PetscFunctionBegin; 769566063dSJacob Faibussowitsch PetscCall(PetscObjectGetComm((PetscObject)pc,&comm)); 77ddc54837SHong Zhang 78ddc54837SHong Zhang /* if pmatrix set by user is sequential then we do not need to gather the parallel matrix */ 799566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_size(comm,&size)); 80ddc54837SHong Zhang if (size == 1) red->useparallelmat = PETSC_FALSE; 811fbd8f88SHong Zhang 824b9ad928SBarry Smith if (!pc->setupcalled) { 831b81debcSHong Zhang PetscInt mloc_sub; 8475024027SHong Zhang if (!red->psubcomm) { /* create red->psubcomm, new ksp and pc over subcomm */ 8575024027SHong Zhang KSP ksp; 869566063dSJacob Faibussowitsch PetscCall(PCRedundantGetKSP(pc,&ksp)); 871b81debcSHong Zhang } 8875024027SHong Zhang subcomm = PetscSubcommChild(red->psubcomm); 891fbd8f88SHong Zhang 901b81debcSHong Zhang if (red->useparallelmat) { 911b81debcSHong Zhang /* grab the parallel matrix and put it into processors of a subcomminicator */ 929566063dSJacob Faibussowitsch PetscCall(MatCreateRedundantMatrix(pc->pmat,red->psubcomm->n,subcomm,MAT_INITIAL_MATRIX,&red->pmats)); 93b85f2e9bSHong Zhang 949566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_size(subcomm,&size)); 95b85f2e9bSHong Zhang if (size > 1) { 9608cecb0aSPierre Jolivet PetscBool foundpack,issbaij; 979566063dSJacob Faibussowitsch PetscCall(PetscObjectTypeCompare((PetscObject)red->pmats,MATMPISBAIJ,&issbaij)); 9808cecb0aSPierre Jolivet if (!issbaij) { 999566063dSJacob Faibussowitsch PetscCall(MatGetFactorAvailable(red->pmats,NULL,MAT_FACTOR_LU,&foundpack)); 10008cecb0aSPierre Jolivet } else { 1019566063dSJacob Faibussowitsch PetscCall(MatGetFactorAvailable(red->pmats,NULL,MAT_FACTOR_CHOLESKY,&foundpack)); 10208cecb0aSPierre Jolivet } 103b85f2e9bSHong Zhang if (!foundpack) { /* reset default ksp and pc */ 1049566063dSJacob Faibussowitsch PetscCall(KSPSetType(red->ksp,KSPGMRES)); 1059566063dSJacob Faibussowitsch PetscCall(PCSetType(red->pc,PCBJACOBI)); 106c1619fb6SBarry Smith } else { 1079566063dSJacob Faibussowitsch PetscCall(PCFactorSetMatSolverType(red->pc,NULL)); 108b85f2e9bSHong Zhang } 109b85f2e9bSHong Zhang } 110b85f2e9bSHong Zhang 1119566063dSJacob Faibussowitsch PetscCall(KSPSetOperators(red->ksp,red->pmats,red->pmats)); 1124b9ad928SBarry Smith 1131b81debcSHong Zhang /* get working vectors xsub and ysub */ 1149566063dSJacob Faibussowitsch PetscCall(MatCreateVecs(red->pmats,&red->xsub,&red->ysub)); 1152fa5cd67SKarl Rupp 1168b96b0d2SHong Zhang /* create working vectors xdup and ydup. 1178b96b0d2SHong Zhang xdup concatenates all xsub's contigously to form a mpi vector over dupcomm (see PetscSubcommCreate_interlaced()) 1188b96b0d2SHong Zhang ydup concatenates all ysub and has empty local arrays because ysub's arrays will be place into it. 119ce94432eSBarry Smith Note: we use communicator dupcomm, not PetscObjectComm((PetscObject)pc)! */ 1209566063dSJacob Faibussowitsch PetscCall(MatGetLocalSize(red->pmats,&mloc_sub,NULL)); 1219566063dSJacob Faibussowitsch PetscCall(VecCreateMPI(PetscSubcommContiguousParent(red->psubcomm),mloc_sub,PETSC_DECIDE,&red->xdup)); 1229566063dSJacob Faibussowitsch PetscCall(VecCreateMPIWithArray(PetscSubcommContiguousParent(red->psubcomm),1,mloc_sub,PETSC_DECIDE,NULL,&red->ydup)); 1233f457be1SHong Zhang 124f68be91cSHong Zhang /* create vecscatters */ 125f68be91cSHong Zhang if (!red->scatterin) { /* efficiency of scatterin is independent from psubcomm_type! */ 1263f457be1SHong Zhang IS is1,is2; 1273f457be1SHong Zhang PetscInt *idx1,*idx2,i,j,k; 12845fc02eaSBarry Smith 1299566063dSJacob Faibussowitsch PetscCall(MatCreateVecs(pc->pmat,&x,NULL)); 1309566063dSJacob Faibussowitsch PetscCall(VecGetSize(x,&M)); 1319566063dSJacob Faibussowitsch PetscCall(VecGetOwnershipRange(x,&mstart,&mend)); 1321b81debcSHong Zhang mlocal = mend - mstart; 1339566063dSJacob Faibussowitsch PetscCall(PetscMalloc2(red->psubcomm->n*mlocal,&idx1,red->psubcomm->n*mlocal,&idx2)); 1343f457be1SHong Zhang j = 0; 1351fbd8f88SHong Zhang for (k=0; k<red->psubcomm->n; k++) { 1363f457be1SHong Zhang for (i=mstart; i<mend; i++) { 1373f457be1SHong Zhang idx1[j] = i; 138ddc54837SHong Zhang idx2[j++] = i + M*k; 1393f457be1SHong Zhang } 1403f457be1SHong Zhang } 1419566063dSJacob Faibussowitsch PetscCall(ISCreateGeneral(comm,red->psubcomm->n*mlocal,idx1,PETSC_COPY_VALUES,&is1)); 1429566063dSJacob Faibussowitsch PetscCall(ISCreateGeneral(comm,red->psubcomm->n*mlocal,idx2,PETSC_COPY_VALUES,&is2)); 1439566063dSJacob Faibussowitsch PetscCall(VecScatterCreate(x,is1,red->xdup,is2,&red->scatterin)); 1449566063dSJacob Faibussowitsch PetscCall(ISDestroy(&is1)); 1459566063dSJacob Faibussowitsch PetscCall(ISDestroy(&is2)); 1463f457be1SHong Zhang 1476909748bSHong Zhang /* Impl below is good for PETSC_SUBCOMM_INTERLACED (no inter-process communication) and PETSC_SUBCOMM_CONTIGUOUS (communication within subcomm) */ 1489566063dSJacob Faibussowitsch PetscCall(ISCreateStride(comm,mlocal,mstart+ red->psubcomm->color*M,1,&is1)); 1499566063dSJacob Faibussowitsch PetscCall(ISCreateStride(comm,mlocal,mstart,1,&is2)); 1509566063dSJacob Faibussowitsch PetscCall(VecScatterCreate(red->xdup,is1,x,is2,&red->scatterout)); 1519566063dSJacob Faibussowitsch PetscCall(ISDestroy(&is1)); 1529566063dSJacob Faibussowitsch PetscCall(ISDestroy(&is2)); 1539566063dSJacob Faibussowitsch PetscCall(PetscFree2(idx1,idx2)); 1549566063dSJacob Faibussowitsch PetscCall(VecDestroy(&x)); 1551b81debcSHong Zhang } 156ab661555SHong Zhang } else { /* !red->useparallelmat */ 1579566063dSJacob Faibussowitsch PetscCall(KSPSetOperators(red->ksp,pc->mat,pc->pmat)); 1581b81debcSHong Zhang } 159ab661555SHong Zhang } else { /* pc->setupcalled */ 1604b9ad928SBarry Smith if (red->useparallelmat) { 161ab661555SHong Zhang MatReuse reuse; 1621b81debcSHong Zhang /* grab the parallel matrix and put it into processors of a subcomminicator */ 1631b81debcSHong Zhang /*--------------------------------------------------------------------------*/ 164ab661555SHong Zhang if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 1654b9ad928SBarry Smith /* destroy old matrices */ 1669566063dSJacob Faibussowitsch PetscCall(MatDestroy(&red->pmats)); 167ab661555SHong Zhang reuse = MAT_INITIAL_MATRIX; 1684b9ad928SBarry Smith } else { 169ab661555SHong Zhang reuse = MAT_REUSE_MATRIX; 170ab661555SHong Zhang } 1719566063dSJacob Faibussowitsch PetscCall(MatCreateRedundantMatrix(pc->pmat,red->psubcomm->n,PetscSubcommChild(red->psubcomm),reuse,&red->pmats)); 1729566063dSJacob Faibussowitsch PetscCall(KSPSetOperators(red->ksp,red->pmats,red->pmats)); 173ab661555SHong Zhang } else { /* !red->useparallelmat */ 1749566063dSJacob Faibussowitsch PetscCall(KSPSetOperators(red->ksp,pc->mat,pc->pmat)); 1754b9ad928SBarry Smith } 176ab661555SHong Zhang } 1771b81debcSHong Zhang 178*1baa6e33SBarry Smith if (pc->setfromoptionscalled) PetscCall(KSPSetFromOptions(red->ksp)); 1799566063dSJacob Faibussowitsch PetscCall(KSPSetUp(red->ksp)); 1804b9ad928SBarry Smith PetscFunctionReturn(0); 1814b9ad928SBarry Smith } 1824b9ad928SBarry Smith 1836849ba73SBarry Smith static PetscErrorCode PCApply_Redundant(PC pc,Vec x,Vec y) 1844b9ad928SBarry Smith { 1854b9ad928SBarry Smith PC_Redundant *red = (PC_Redundant*)pc->data; 1863f457be1SHong Zhang PetscScalar *array; 1874b9ad928SBarry Smith 1884b9ad928SBarry Smith PetscFunctionBegin; 189ddc54837SHong Zhang if (!red->useparallelmat) { 1909566063dSJacob Faibussowitsch PetscCall(KSPSolve(red->ksp,x,y)); 1919566063dSJacob Faibussowitsch PetscCall(KSPCheckSolve(red->ksp,pc,y)); 192ddc54837SHong Zhang PetscFunctionReturn(0); 193ddc54837SHong Zhang } 194ddc54837SHong Zhang 1953f457be1SHong Zhang /* scatter x to xdup */ 1969566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD)); 1979566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD)); 1983f457be1SHong Zhang 1993f457be1SHong Zhang /* place xdup's local array into xsub */ 2009566063dSJacob Faibussowitsch PetscCall(VecGetArray(red->xdup,&array)); 2019566063dSJacob Faibussowitsch PetscCall(VecPlaceArray(red->xsub,(const PetscScalar*)array)); 2024b9ad928SBarry Smith 2034b9ad928SBarry Smith /* apply preconditioner on each processor */ 2049566063dSJacob Faibussowitsch PetscCall(KSPSolve(red->ksp,red->xsub,red->ysub)); 2059566063dSJacob Faibussowitsch PetscCall(KSPCheckSolve(red->ksp,pc,red->ysub)); 2069566063dSJacob Faibussowitsch PetscCall(VecResetArray(red->xsub)); 2079566063dSJacob Faibussowitsch PetscCall(VecRestoreArray(red->xdup,&array)); 2084b9ad928SBarry Smith 2093f457be1SHong Zhang /* place ysub's local array into ydup */ 2109566063dSJacob Faibussowitsch PetscCall(VecGetArray(red->ysub,&array)); 2119566063dSJacob Faibussowitsch PetscCall(VecPlaceArray(red->ydup,(const PetscScalar*)array)); 2123f457be1SHong Zhang 2133f457be1SHong Zhang /* scatter ydup to y */ 2149566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD)); 2159566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD)); 2169566063dSJacob Faibussowitsch PetscCall(VecResetArray(red->ydup)); 2179566063dSJacob Faibussowitsch PetscCall(VecRestoreArray(red->ysub,&array)); 2184b9ad928SBarry Smith PetscFunctionReturn(0); 2194b9ad928SBarry Smith } 2204b9ad928SBarry Smith 221d88bfacbSStefano Zampini static PetscErrorCode PCApplyTranspose_Redundant(PC pc,Vec x,Vec y) 222d88bfacbSStefano Zampini { 223d88bfacbSStefano Zampini PC_Redundant *red = (PC_Redundant*)pc->data; 224d88bfacbSStefano Zampini PetscScalar *array; 225d88bfacbSStefano Zampini 226d88bfacbSStefano Zampini PetscFunctionBegin; 227d88bfacbSStefano Zampini if (!red->useparallelmat) { 2289566063dSJacob Faibussowitsch PetscCall(KSPSolveTranspose(red->ksp,x,y)); 2299566063dSJacob Faibussowitsch PetscCall(KSPCheckSolve(red->ksp,pc,y)); 230d88bfacbSStefano Zampini PetscFunctionReturn(0); 231d88bfacbSStefano Zampini } 232d88bfacbSStefano Zampini 233d88bfacbSStefano Zampini /* scatter x to xdup */ 2349566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD)); 2359566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD)); 236d88bfacbSStefano Zampini 237d88bfacbSStefano Zampini /* place xdup's local array into xsub */ 2389566063dSJacob Faibussowitsch PetscCall(VecGetArray(red->xdup,&array)); 2399566063dSJacob Faibussowitsch PetscCall(VecPlaceArray(red->xsub,(const PetscScalar*)array)); 240d88bfacbSStefano Zampini 241d88bfacbSStefano Zampini /* apply preconditioner on each processor */ 2429566063dSJacob Faibussowitsch PetscCall(KSPSolveTranspose(red->ksp,red->xsub,red->ysub)); 2439566063dSJacob Faibussowitsch PetscCall(KSPCheckSolve(red->ksp,pc,red->ysub)); 2449566063dSJacob Faibussowitsch PetscCall(VecResetArray(red->xsub)); 2459566063dSJacob Faibussowitsch PetscCall(VecRestoreArray(red->xdup,&array)); 246d88bfacbSStefano Zampini 247d88bfacbSStefano Zampini /* place ysub's local array into ydup */ 2489566063dSJacob Faibussowitsch PetscCall(VecGetArray(red->ysub,&array)); 2499566063dSJacob Faibussowitsch PetscCall(VecPlaceArray(red->ydup,(const PetscScalar*)array)); 250d88bfacbSStefano Zampini 251d88bfacbSStefano Zampini /* scatter ydup to y */ 2529566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD)); 2539566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD)); 2549566063dSJacob Faibussowitsch PetscCall(VecResetArray(red->ydup)); 2559566063dSJacob Faibussowitsch PetscCall(VecRestoreArray(red->ysub,&array)); 256d88bfacbSStefano Zampini PetscFunctionReturn(0); 257d88bfacbSStefano Zampini } 258d88bfacbSStefano Zampini 2591ea5a559SBarry Smith static PetscErrorCode PCReset_Redundant(PC pc) 2604b9ad928SBarry Smith { 2614b9ad928SBarry Smith PC_Redundant *red = (PC_Redundant*)pc->data; 2624b9ad928SBarry Smith 2634b9ad928SBarry Smith PetscFunctionBegin; 2641b81debcSHong Zhang if (red->useparallelmat) { 2659566063dSJacob Faibussowitsch PetscCall(VecScatterDestroy(&red->scatterin)); 2669566063dSJacob Faibussowitsch PetscCall(VecScatterDestroy(&red->scatterout)); 2679566063dSJacob Faibussowitsch PetscCall(VecDestroy(&red->ysub)); 2689566063dSJacob Faibussowitsch PetscCall(VecDestroy(&red->xsub)); 2699566063dSJacob Faibussowitsch PetscCall(VecDestroy(&red->xdup)); 2709566063dSJacob Faibussowitsch PetscCall(VecDestroy(&red->ydup)); 2711b81debcSHong Zhang } 2729566063dSJacob Faibussowitsch PetscCall(MatDestroy(&red->pmats)); 2739566063dSJacob Faibussowitsch PetscCall(KSPReset(red->ksp)); 2741ea5a559SBarry Smith PetscFunctionReturn(0); 2751ea5a559SBarry Smith } 2761ea5a559SBarry Smith 2771ea5a559SBarry Smith static PetscErrorCode PCDestroy_Redundant(PC pc) 2781ea5a559SBarry Smith { 2791ea5a559SBarry Smith PC_Redundant *red = (PC_Redundant*)pc->data; 2801ea5a559SBarry Smith 2811ea5a559SBarry Smith PetscFunctionBegin; 2829566063dSJacob Faibussowitsch PetscCall(PCReset_Redundant(pc)); 2839566063dSJacob Faibussowitsch PetscCall(KSPDestroy(&red->ksp)); 2849566063dSJacob Faibussowitsch PetscCall(PetscSubcommDestroy(&red->psubcomm)); 2852e956fe4SStefano Zampini PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCRedundantSetScatter_C",NULL)); 2862e956fe4SStefano Zampini PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCRedundantSetNumber_C",NULL)); 2872e956fe4SStefano Zampini PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCRedundantGetKSP_C",NULL)); 2882e956fe4SStefano Zampini PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCRedundantGetOperators_C",NULL)); 2892e956fe4SStefano Zampini PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftType_C",NULL)); 2909566063dSJacob Faibussowitsch PetscCall(PetscFree(pc->data)); 2914b9ad928SBarry Smith PetscFunctionReturn(0); 2924b9ad928SBarry Smith } 2934b9ad928SBarry Smith 2944416b707SBarry Smith static PetscErrorCode PCSetFromOptions_Redundant(PetscOptionItems *PetscOptionsObject,PC pc) 2954b9ad928SBarry Smith { 296a98ce0f4SHong Zhang PC_Redundant *red = (PC_Redundant*)pc->data; 297a98ce0f4SHong Zhang 2984b9ad928SBarry Smith PetscFunctionBegin; 299d0609cedSBarry Smith PetscOptionsHeadBegin(PetscOptionsObject,"Redundant options"); 3009566063dSJacob Faibussowitsch PetscCall(PetscOptionsInt("-pc_redundant_number","Number of redundant pc","PCRedundantSetNumber",red->nsubcomm,&red->nsubcomm,NULL)); 301d0609cedSBarry Smith PetscOptionsHeadEnd(); 3024b9ad928SBarry Smith PetscFunctionReturn(0); 3034b9ad928SBarry Smith } 3044b9ad928SBarry Smith 305f7a08781SBarry Smith static PetscErrorCode PCRedundantSetNumber_Redundant(PC pc,PetscInt nreds) 30609a6bc64SHong Zhang { 30709a6bc64SHong Zhang PC_Redundant *red = (PC_Redundant*)pc->data; 30809a6bc64SHong Zhang 30909a6bc64SHong Zhang PetscFunctionBegin; 31009a6bc64SHong Zhang red->nsubcomm = nreds; 31109a6bc64SHong Zhang PetscFunctionReturn(0); 31209a6bc64SHong Zhang } 31309a6bc64SHong Zhang 31409a6bc64SHong Zhang /*@ 31509a6bc64SHong Zhang PCRedundantSetNumber - Sets the number of redundant preconditioner contexts. 31609a6bc64SHong Zhang 3173f9fe445SBarry Smith Logically Collective on PC 31809a6bc64SHong Zhang 31909a6bc64SHong Zhang Input Parameters: 32009a6bc64SHong Zhang + pc - the preconditioner context 3219b21d695SBarry Smith - nredundant - number of redundant preconditioner contexts; for example if you are using 64 MPI processes and 3229b21d695SBarry Smith use an nredundant of 4 there will be 4 parallel solves each on 16 = 64/4 processes. 32309a6bc64SHong Zhang 32409a6bc64SHong Zhang Level: advanced 32509a6bc64SHong Zhang 32609a6bc64SHong Zhang @*/ 3277087cfbeSBarry Smith PetscErrorCode PCRedundantSetNumber(PC pc,PetscInt nredundant) 32809a6bc64SHong Zhang { 32909a6bc64SHong Zhang PetscFunctionBegin; 3300700a824SBarry Smith PetscValidHeaderSpecific(pc,PC_CLASSID,1); 33163a3b9bcSJacob Faibussowitsch PetscCheck(nredundant > 0,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "num of redundant pc %" PetscInt_FMT " must be positive",nredundant); 332cac4c232SBarry Smith PetscTryMethod(pc,"PCRedundantSetNumber_C",(PC,PetscInt),(pc,nredundant)); 33309a6bc64SHong Zhang PetscFunctionReturn(0); 33409a6bc64SHong Zhang } 33509a6bc64SHong Zhang 336f7a08781SBarry Smith static PetscErrorCode PCRedundantSetScatter_Redundant(PC pc,VecScatter in,VecScatter out) 3374b9ad928SBarry Smith { 3384b9ad928SBarry Smith PC_Redundant *red = (PC_Redundant*)pc->data; 3394b9ad928SBarry Smith 3404b9ad928SBarry Smith PetscFunctionBegin; 3419566063dSJacob Faibussowitsch PetscCall(PetscObjectReference((PetscObject)in)); 3429566063dSJacob Faibussowitsch PetscCall(VecScatterDestroy(&red->scatterin)); 3432fa5cd67SKarl Rupp 344c3122656SLisandro Dalcin red->scatterin = in; 3452fa5cd67SKarl Rupp 3469566063dSJacob Faibussowitsch PetscCall(PetscObjectReference((PetscObject)out)); 3479566063dSJacob Faibussowitsch PetscCall(VecScatterDestroy(&red->scatterout)); 348c3122656SLisandro Dalcin red->scatterout = out; 3494b9ad928SBarry Smith PetscFunctionReturn(0); 3504b9ad928SBarry Smith } 3514b9ad928SBarry Smith 3524b9ad928SBarry Smith /*@ 3534b9ad928SBarry Smith PCRedundantSetScatter - Sets the scatter used to copy values into the 3544b9ad928SBarry Smith redundant local solve and the scatter to move them back into the global 3554b9ad928SBarry Smith vector. 3564b9ad928SBarry Smith 3573f9fe445SBarry Smith Logically Collective on PC 3584b9ad928SBarry Smith 3594b9ad928SBarry Smith Input Parameters: 3604b9ad928SBarry Smith + pc - the preconditioner context 3614b9ad928SBarry Smith . in - the scatter to move the values in 3624b9ad928SBarry Smith - out - the scatter to move them out 3634b9ad928SBarry Smith 3644b9ad928SBarry Smith Level: advanced 3654b9ad928SBarry Smith 3664b9ad928SBarry Smith @*/ 3677087cfbeSBarry Smith PetscErrorCode PCRedundantSetScatter(PC pc,VecScatter in,VecScatter out) 3684b9ad928SBarry Smith { 3694b9ad928SBarry Smith PetscFunctionBegin; 3700700a824SBarry Smith PetscValidHeaderSpecific(pc,PC_CLASSID,1); 37197929ea7SJunchao Zhang PetscValidHeaderSpecific(in,PETSCSF_CLASSID,2); 37297929ea7SJunchao Zhang PetscValidHeaderSpecific(out,PETSCSF_CLASSID,3); 373cac4c232SBarry Smith PetscTryMethod(pc,"PCRedundantSetScatter_C",(PC,VecScatter,VecScatter),(pc,in,out)); 3744b9ad928SBarry Smith PetscFunctionReturn(0); 3754b9ad928SBarry Smith } 3764b9ad928SBarry Smith 377f7a08781SBarry Smith static PetscErrorCode PCRedundantGetKSP_Redundant(PC pc,KSP *innerksp) 3784b9ad928SBarry Smith { 3794b9ad928SBarry Smith PC_Redundant *red = (PC_Redundant*)pc->data; 38075024027SHong Zhang MPI_Comm comm,subcomm; 38175024027SHong Zhang const char *prefix; 38208cecb0aSPierre Jolivet PetscBool issbaij; 3834b9ad928SBarry Smith 3844b9ad928SBarry Smith PetscFunctionBegin; 38575024027SHong Zhang if (!red->psubcomm) { 3869566063dSJacob Faibussowitsch PetscCall(PCGetOptionsPrefix(pc,&prefix)); 387e5acf8a4SHong Zhang 3889566063dSJacob Faibussowitsch PetscCall(PetscObjectGetComm((PetscObject)pc,&comm)); 3899566063dSJacob Faibussowitsch PetscCall(PetscSubcommCreate(comm,&red->psubcomm)); 3909566063dSJacob Faibussowitsch PetscCall(PetscSubcommSetNumber(red->psubcomm,red->nsubcomm)); 3919566063dSJacob Faibussowitsch PetscCall(PetscSubcommSetType(red->psubcomm,PETSC_SUBCOMM_CONTIGUOUS)); 392e5acf8a4SHong Zhang 3939566063dSJacob Faibussowitsch PetscCall(PetscSubcommSetOptionsPrefix(red->psubcomm,prefix)); 3949566063dSJacob Faibussowitsch PetscCall(PetscSubcommSetFromOptions(red->psubcomm)); 3959566063dSJacob Faibussowitsch PetscCall(PetscLogObjectMemory((PetscObject)pc,sizeof(PetscSubcomm))); 39675024027SHong Zhang 39775024027SHong Zhang /* create a new PC that processors in each subcomm have copy of */ 39875024027SHong Zhang subcomm = PetscSubcommChild(red->psubcomm); 39975024027SHong Zhang 4009566063dSJacob Faibussowitsch PetscCall(KSPCreate(subcomm,&red->ksp)); 4019566063dSJacob Faibussowitsch PetscCall(KSPSetErrorIfNotConverged(red->ksp,pc->erroriffailure)); 4029566063dSJacob Faibussowitsch PetscCall(PetscObjectIncrementTabLevel((PetscObject)red->ksp,(PetscObject)pc,1)); 4039566063dSJacob Faibussowitsch PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)red->ksp)); 4049566063dSJacob Faibussowitsch PetscCall(KSPSetType(red->ksp,KSPPREONLY)); 4059566063dSJacob Faibussowitsch PetscCall(KSPGetPC(red->ksp,&red->pc)); 4069566063dSJacob Faibussowitsch PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat,MATSEQSBAIJ,&issbaij)); 40708cecb0aSPierre Jolivet if (!issbaij) { 4089566063dSJacob Faibussowitsch PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat,MATMPISBAIJ,&issbaij)); 40908cecb0aSPierre Jolivet } 41008cecb0aSPierre Jolivet if (!issbaij) { 4119566063dSJacob Faibussowitsch PetscCall(PCSetType(red->pc,PCLU)); 41208cecb0aSPierre Jolivet } else { 4139566063dSJacob Faibussowitsch PetscCall(PCSetType(red->pc,PCCHOLESKY)); 41408cecb0aSPierre Jolivet } 415753b7fb9SBarry Smith if (red->shifttypeset) { 4169566063dSJacob Faibussowitsch PetscCall(PCFactorSetShiftType(red->pc,red->shifttype)); 417753b7fb9SBarry Smith red->shifttypeset = PETSC_FALSE; 418753b7fb9SBarry Smith } 4199566063dSJacob Faibussowitsch PetscCall(KSPSetOptionsPrefix(red->ksp,prefix)); 4209566063dSJacob Faibussowitsch PetscCall(KSPAppendOptionsPrefix(red->ksp,"redundant_")); 42175024027SHong Zhang } 42283ab6a24SBarry Smith *innerksp = red->ksp; 4234b9ad928SBarry Smith PetscFunctionReturn(0); 4244b9ad928SBarry Smith } 4254b9ad928SBarry Smith 4264b9ad928SBarry Smith /*@ 42783ab6a24SBarry Smith PCRedundantGetKSP - Gets the less parallel KSP created by the redundant PC. 4284b9ad928SBarry Smith 4294b9ad928SBarry Smith Not Collective 4304b9ad928SBarry Smith 4314b9ad928SBarry Smith Input Parameter: 4324b9ad928SBarry Smith . pc - the preconditioner context 4334b9ad928SBarry Smith 4344b9ad928SBarry Smith Output Parameter: 43583ab6a24SBarry Smith . innerksp - the KSP on the smaller set of processes 4364b9ad928SBarry Smith 4374b9ad928SBarry Smith Level: advanced 4384b9ad928SBarry Smith 4394b9ad928SBarry Smith @*/ 44083ab6a24SBarry Smith PetscErrorCode PCRedundantGetKSP(PC pc,KSP *innerksp) 4414b9ad928SBarry Smith { 4424b9ad928SBarry Smith PetscFunctionBegin; 4430700a824SBarry Smith PetscValidHeaderSpecific(pc,PC_CLASSID,1); 44483ab6a24SBarry Smith PetscValidPointer(innerksp,2); 445cac4c232SBarry Smith PetscUseMethod(pc,"PCRedundantGetKSP_C",(PC,KSP*),(pc,innerksp)); 4464b9ad928SBarry Smith PetscFunctionReturn(0); 4474b9ad928SBarry Smith } 4484b9ad928SBarry Smith 449f7a08781SBarry Smith static PetscErrorCode PCRedundantGetOperators_Redundant(PC pc,Mat *mat,Mat *pmat) 4504b9ad928SBarry Smith { 4514b9ad928SBarry Smith PC_Redundant *red = (PC_Redundant*)pc->data; 4524b9ad928SBarry Smith 4534b9ad928SBarry Smith PetscFunctionBegin; 454b3804887SHong Zhang if (mat) *mat = red->pmats; 455b3804887SHong Zhang if (pmat) *pmat = red->pmats; 4564b9ad928SBarry Smith PetscFunctionReturn(0); 4574b9ad928SBarry Smith } 4584b9ad928SBarry Smith 4594b9ad928SBarry Smith /*@ 4604b9ad928SBarry Smith PCRedundantGetOperators - gets the sequential matrix and preconditioner matrix 4614b9ad928SBarry Smith 4624b9ad928SBarry Smith Not Collective 4634b9ad928SBarry Smith 4644b9ad928SBarry Smith Input Parameter: 4654b9ad928SBarry Smith . pc - the preconditioner context 4664b9ad928SBarry Smith 4674b9ad928SBarry Smith Output Parameters: 4684b9ad928SBarry Smith + mat - the matrix 4694b9ad928SBarry Smith - pmat - the (possibly different) preconditioner matrix 4704b9ad928SBarry Smith 4714b9ad928SBarry Smith Level: advanced 4724b9ad928SBarry Smith 4734b9ad928SBarry Smith @*/ 4747087cfbeSBarry Smith PetscErrorCode PCRedundantGetOperators(PC pc,Mat *mat,Mat *pmat) 4754b9ad928SBarry Smith { 4764b9ad928SBarry Smith PetscFunctionBegin; 4770700a824SBarry Smith PetscValidHeaderSpecific(pc,PC_CLASSID,1); 4784482741eSBarry Smith if (mat) PetscValidPointer(mat,2); 4794482741eSBarry Smith if (pmat) PetscValidPointer(pmat,3); 480cac4c232SBarry Smith PetscUseMethod(pc,"PCRedundantGetOperators_C",(PC,Mat*,Mat*),(pc,mat,pmat)); 4814b9ad928SBarry Smith PetscFunctionReturn(0); 4824b9ad928SBarry Smith } 4834b9ad928SBarry Smith 4844b9ad928SBarry Smith /* -------------------------------------------------------------------------------------*/ 48537a17b4dSBarry Smith /*MC 48683ab6a24SBarry Smith PCREDUNDANT - Runs a KSP solver with preconditioner for the entire problem on subgroups of processors 48737a17b4dSBarry Smith 48883ab6a24SBarry Smith Options for the redundant preconditioners can be set with -redundant_pc_xxx for the redundant KSP with -redundant_ksp_xxx 48937a17b4dSBarry Smith 49009391456SBarry Smith Options Database: 4919b21d695SBarry Smith . -pc_redundant_number <n> - number of redundant solves, for example if you are using 64 MPI processes and 4929b21d695SBarry Smith use an n of 4 there will be 4 parallel solves each on 16 = 64/4 processes. 49309391456SBarry Smith 49437a17b4dSBarry Smith Level: intermediate 49537a17b4dSBarry Smith 49695452b02SPatrick Sanan Notes: 49708cecb0aSPierre Jolivet The default KSP is preonly and the default PC is LU or CHOLESKY if Pmat is of type MATSBAIJ. 49883ab6a24SBarry Smith 499753b7fb9SBarry Smith PCFactorSetShiftType() applied to this PC will convey they shift type into the inner PC if it is factorization based. 500753b7fb9SBarry Smith 50195452b02SPatrick Sanan Developer Notes: 50295452b02SPatrick Sanan Note that PCSetInitialGuessNonzero() is not used by this class but likely should be. 5039cfaa89bSBarry Smith 504db781477SPatrick Sanan .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PCRedundantSetScatter()`, 505db781477SPatrick Sanan `PCRedundantGetKSP()`, `PCRedundantGetOperators()`, `PCRedundantSetNumber()` 50637a17b4dSBarry Smith M*/ 50737a17b4dSBarry Smith 5088cc058d9SJed Brown PETSC_EXTERN PetscErrorCode PCCreate_Redundant(PC pc) 5094b9ad928SBarry Smith { 5104b9ad928SBarry Smith PC_Redundant *red; 51169db28dcSHong Zhang PetscMPIInt size; 5123f457be1SHong Zhang 5134b9ad928SBarry Smith PetscFunctionBegin; 5149566063dSJacob Faibussowitsch PetscCall(PetscNewLog(pc,&red)); 5159566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size)); 5162fa5cd67SKarl Rupp 51769db28dcSHong Zhang red->nsubcomm = size; 5184b9ad928SBarry Smith red->useparallelmat = PETSC_TRUE; 5191fbd8f88SHong Zhang pc->data = (void*)red; 5204b9ad928SBarry Smith 5214b9ad928SBarry Smith pc->ops->apply = PCApply_Redundant; 522d88bfacbSStefano Zampini pc->ops->applytranspose = PCApplyTranspose_Redundant; 5234b9ad928SBarry Smith pc->ops->setup = PCSetUp_Redundant; 5244b9ad928SBarry Smith pc->ops->destroy = PCDestroy_Redundant; 5251ea5a559SBarry Smith pc->ops->reset = PCReset_Redundant; 5264b9ad928SBarry Smith pc->ops->setfromoptions = PCSetFromOptions_Redundant; 5274b9ad928SBarry Smith pc->ops->view = PCView_Redundant; 5282fa5cd67SKarl Rupp 5299566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCRedundantSetScatter_C",PCRedundantSetScatter_Redundant)); 5309566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCRedundantSetNumber_C",PCRedundantSetNumber_Redundant)); 5319566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCRedundantGetKSP_C",PCRedundantGetKSP_Redundant)); 5329566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCRedundantGetOperators_C",PCRedundantGetOperators_Redundant)); 5339566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftType_C",PCFactorSetShiftType_Redundant)); 5344b9ad928SBarry Smith PetscFunctionReturn(0); 5354b9ad928SBarry Smith } 536