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 PetscErrorCode ierr; 26753b7fb9SBarry Smith 27753b7fb9SBarry Smith PetscFunctionBegin; 28753b7fb9SBarry Smith if (red->ksp) { 29753b7fb9SBarry Smith PC pc; 30753b7fb9SBarry Smith ierr = KSPGetPC(red->ksp,&pc);CHKERRQ(ierr); 31753b7fb9SBarry Smith ierr = PCFactorSetShiftType(pc,shifttype);CHKERRQ(ierr); 32753b7fb9SBarry Smith } else { 33753b7fb9SBarry Smith red->shifttypeset = PETSC_TRUE; 34753b7fb9SBarry Smith red->shifttype = shifttype; 35753b7fb9SBarry Smith } 36753b7fb9SBarry Smith PetscFunctionReturn(0); 37753b7fb9SBarry Smith } 38753b7fb9SBarry Smith 396849ba73SBarry Smith static PetscErrorCode PCView_Redundant(PC pc,PetscViewer viewer) 404b9ad928SBarry Smith { 414b9ad928SBarry Smith PC_Redundant *red = (PC_Redundant*)pc->data; 42dfbe8321SBarry Smith PetscErrorCode ierr; 43ace3abfcSBarry Smith PetscBool iascii,isstring; 4403ccd0b4SBarry Smith PetscViewer subviewer; 454b9ad928SBarry Smith 464b9ad928SBarry Smith PetscFunctionBegin; 47251f4c67SDmitry Karpeev ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 48251f4c67SDmitry Karpeev ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);CHKERRQ(ierr); 4932077d6dSBarry Smith if (iascii) { 5003ccd0b4SBarry Smith if (!red->psubcomm) { 51*efd4aadfSBarry Smith ierr = PetscViewerASCIIPrintf(viewer," Not yet setup\n");CHKERRQ(ierr); 5203ccd0b4SBarry Smith } else { 53*efd4aadfSBarry Smith ierr = PetscViewerASCIIPrintf(viewer," First (color=0) of %D PCs follows\n",red->nsubcomm);CHKERRQ(ierr); 543f08860eSBarry Smith ierr = PetscViewerGetSubViewer(viewer,((PetscObject)red->pc)->comm,&subviewer);CHKERRQ(ierr); 55f5dd71faSBarry Smith if (!red->psubcomm->color) { /* only view first redundant pc */ 561575c14dSBarry Smith ierr = PetscViewerASCIIPushTab(subviewer);CHKERRQ(ierr); 573e065800SHong Zhang ierr = KSPView(red->ksp,subviewer);CHKERRQ(ierr); 581575c14dSBarry Smith ierr = PetscViewerASCIIPopTab(subviewer);CHKERRQ(ierr); 594b9ad928SBarry Smith } 603f08860eSBarry Smith ierr = PetscViewerRestoreSubViewer(viewer,((PetscObject)red->pc)->comm,&subviewer);CHKERRQ(ierr); 614b9ad928SBarry Smith } 6203ccd0b4SBarry Smith } else if (isstring) { 6303ccd0b4SBarry Smith ierr = PetscViewerStringSPrintf(viewer," Redundant solver preconditioner");CHKERRQ(ierr); 644b9ad928SBarry Smith } 654b9ad928SBarry Smith PetscFunctionReturn(0); 664b9ad928SBarry Smith } 674b9ad928SBarry Smith 686849ba73SBarry Smith static PetscErrorCode PCSetUp_Redundant(PC pc) 694b9ad928SBarry Smith { 704b9ad928SBarry Smith PC_Redundant *red = (PC_Redundant*)pc->data; 71dfbe8321SBarry Smith PetscErrorCode ierr; 721b81debcSHong Zhang PetscInt mstart,mend,mlocal,M; 7313f74950SBarry Smith PetscMPIInt size; 74ce94432eSBarry Smith MPI_Comm comm,subcomm; 75ddc54837SHong Zhang Vec x; 763f457be1SHong Zhang 774b9ad928SBarry Smith PetscFunctionBegin; 78ce94432eSBarry Smith ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 79ddc54837SHong Zhang 80ddc54837SHong Zhang /* if pmatrix set by user is sequential then we do not need to gather the parallel matrix */ 81ddc54837SHong Zhang ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 82ddc54837SHong Zhang if (size == 1) red->useparallelmat = PETSC_FALSE; 831fbd8f88SHong Zhang 844b9ad928SBarry Smith if (!pc->setupcalled) { 851b81debcSHong Zhang PetscInt mloc_sub; 8675024027SHong Zhang if (!red->psubcomm) { /* create red->psubcomm, new ksp and pc over subcomm */ 8775024027SHong Zhang KSP ksp; 8875024027SHong Zhang ierr = PCRedundantGetKSP(pc,&ksp);CHKERRQ(ierr); 891b81debcSHong Zhang } 9075024027SHong Zhang subcomm = PetscSubcommChild(red->psubcomm); 911fbd8f88SHong Zhang 921b81debcSHong Zhang if (red->useparallelmat) { 931b81debcSHong Zhang /* grab the parallel matrix and put it into processors of a subcomminicator */ 9453cd1579SHong Zhang ierr = MatCreateRedundantMatrix(pc->pmat,red->psubcomm->n,subcomm,MAT_INITIAL_MATRIX,&red->pmats);CHKERRQ(ierr); 95b85f2e9bSHong Zhang 96b85f2e9bSHong Zhang ierr = MPI_Comm_size(subcomm,&size);CHKERRQ(ierr); 97b85f2e9bSHong Zhang if (size > 1) { 98b85f2e9bSHong Zhang PetscBool foundpack; 99b85f2e9bSHong Zhang ierr = MatGetFactorAvailable(red->pmats,NULL,MAT_FACTOR_LU,&foundpack);CHKERRQ(ierr); 100b85f2e9bSHong Zhang if (!foundpack) { /* reset default ksp and pc */ 101b85f2e9bSHong Zhang ierr = KSPSetType(red->ksp,KSPGMRES);CHKERRQ(ierr); 102b85f2e9bSHong Zhang ierr = PCSetType(red->pc,PCBJACOBI);CHKERRQ(ierr); 103c1619fb6SBarry Smith } else { 104c1619fb6SBarry Smith ierr = PCFactorSetMatSolverPackage(red->pc,NULL);CHKERRQ(ierr); 105b85f2e9bSHong Zhang } 106b85f2e9bSHong Zhang } 107b85f2e9bSHong Zhang 10823ee1639SBarry Smith ierr = KSPSetOperators(red->ksp,red->pmats,red->pmats);CHKERRQ(ierr); 1094b9ad928SBarry Smith 1101b81debcSHong Zhang /* get working vectors xsub and ysub */ 1112a7a6963SBarry Smith ierr = MatCreateVecs(red->pmats,&red->xsub,&red->ysub);CHKERRQ(ierr); 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)! */ 1171b81debcSHong Zhang ierr = MatGetLocalSize(red->pmats,&mloc_sub,NULL);CHKERRQ(ierr); 11836be1a5eSBarry Smith ierr = VecCreateMPI(PetscSubcommContiguousParent(red->psubcomm),mloc_sub,PETSC_DECIDE,&red->xdup);CHKERRQ(ierr); 11936be1a5eSBarry Smith ierr = VecCreateMPIWithArray(PetscSubcommContiguousParent(red->psubcomm),1,mloc_sub,PETSC_DECIDE,NULL,&red->ydup);CHKERRQ(ierr); 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 1262a7a6963SBarry Smith ierr = MatCreateVecs(pc->pmat,&x,0);CHKERRQ(ierr); 1271b81debcSHong Zhang ierr = VecGetSize(x,&M);CHKERRQ(ierr); 1281b81debcSHong Zhang ierr = VecGetOwnershipRange(x,&mstart,&mend);CHKERRQ(ierr); 1291b81debcSHong Zhang mlocal = mend - mstart; 130dcca6d9dSJed Brown ierr = PetscMalloc2(red->psubcomm->n*mlocal,&idx1,red->psubcomm->n*mlocal,&idx2);CHKERRQ(ierr); 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 } 13870b3c8c7SBarry Smith ierr = ISCreateGeneral(comm,red->psubcomm->n*mlocal,idx1,PETSC_COPY_VALUES,&is1);CHKERRQ(ierr); 13970b3c8c7SBarry Smith ierr = ISCreateGeneral(comm,red->psubcomm->n*mlocal,idx2,PETSC_COPY_VALUES,&is2);CHKERRQ(ierr); 140ddc54837SHong Zhang ierr = VecScatterCreate(x,is1,red->xdup,is2,&red->scatterin);CHKERRQ(ierr); 141fcfd50ebSBarry Smith ierr = ISDestroy(&is1);CHKERRQ(ierr); 142fcfd50ebSBarry Smith ierr = ISDestroy(&is2);CHKERRQ(ierr); 1433f457be1SHong Zhang 1446909748bSHong Zhang /* Impl below is good for PETSC_SUBCOMM_INTERLACED (no inter-process communication) and PETSC_SUBCOMM_CONTIGUOUS (communication within subcomm) */ 145ddc54837SHong Zhang ierr = ISCreateStride(comm,mlocal,mstart+ red->psubcomm->color*M,1,&is1);CHKERRQ(ierr); 1463f457be1SHong Zhang ierr = ISCreateStride(comm,mlocal,mstart,1,&is2);CHKERRQ(ierr); 147ddc54837SHong Zhang ierr = VecScatterCreate(red->xdup,is1,x,is2,&red->scatterout);CHKERRQ(ierr); 148fcfd50ebSBarry Smith ierr = ISDestroy(&is1);CHKERRQ(ierr); 149fcfd50ebSBarry Smith ierr = ISDestroy(&is2);CHKERRQ(ierr); 1501d79065fSBarry Smith ierr = PetscFree2(idx1,idx2);CHKERRQ(ierr); 151ddc54837SHong Zhang ierr = VecDestroy(&x);CHKERRQ(ierr); 1521b81debcSHong Zhang } 153ab661555SHong Zhang } else { /* !red->useparallelmat */ 15423ee1639SBarry Smith ierr = KSPSetOperators(red->ksp,pc->mat,pc->pmat);CHKERRQ(ierr); 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 */ 1636bf464f9SBarry Smith ierr = MatDestroy(&red->pmats);CHKERRQ(ierr); 164ab661555SHong Zhang reuse = MAT_INITIAL_MATRIX; 1654b9ad928SBarry Smith } else { 166ab661555SHong Zhang reuse = MAT_REUSE_MATRIX; 167ab661555SHong Zhang } 168306c2d5bSBarry Smith ierr = MatCreateRedundantMatrix(pc->pmat,red->psubcomm->n,PetscSubcommChild(red->psubcomm),reuse,&red->pmats);CHKERRQ(ierr); 16923ee1639SBarry Smith ierr = KSPSetOperators(red->ksp,red->pmats,red->pmats);CHKERRQ(ierr); 170ab661555SHong Zhang } else { /* !red->useparallelmat */ 17123ee1639SBarry Smith ierr = KSPSetOperators(red->ksp,pc->mat,pc->pmat);CHKERRQ(ierr); 1724b9ad928SBarry Smith } 173ab661555SHong Zhang } 1741b81debcSHong Zhang 1750c24e6a1SHong Zhang if (pc->setfromoptionscalled) { 1763e065800SHong Zhang ierr = KSPSetFromOptions(red->ksp);CHKERRQ(ierr); 1770c24e6a1SHong Zhang } 1783e065800SHong Zhang ierr = KSPSetUp(red->ksp);CHKERRQ(ierr); 1794b9ad928SBarry Smith PetscFunctionReturn(0); 1804b9ad928SBarry Smith } 1814b9ad928SBarry Smith 1826849ba73SBarry Smith static PetscErrorCode PCApply_Redundant(PC pc,Vec x,Vec y) 1834b9ad928SBarry Smith { 1844b9ad928SBarry Smith PC_Redundant *red = (PC_Redundant*)pc->data; 185dfbe8321SBarry Smith PetscErrorCode ierr; 1863f457be1SHong Zhang PetscScalar *array; 1874b9ad928SBarry Smith 1884b9ad928SBarry Smith PetscFunctionBegin; 189ddc54837SHong Zhang if (!red->useparallelmat) { 190ddc54837SHong Zhang ierr = KSPSolve(red->ksp,x,y);CHKERRQ(ierr); 191ddc54837SHong Zhang PetscFunctionReturn(0); 192ddc54837SHong Zhang } 193ddc54837SHong Zhang 1943f457be1SHong Zhang /* scatter x to xdup */ 195ca9f406cSSatish Balay ierr = VecScatterBegin(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 196ca9f406cSSatish Balay ierr = VecScatterEnd(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1973f457be1SHong Zhang 1983f457be1SHong Zhang /* place xdup's local array into xsub */ 1993f457be1SHong Zhang ierr = VecGetArray(red->xdup,&array);CHKERRQ(ierr); 2003f457be1SHong Zhang ierr = VecPlaceArray(red->xsub,(const PetscScalar*)array);CHKERRQ(ierr); 2014b9ad928SBarry Smith 2024b9ad928SBarry Smith /* apply preconditioner on each processor */ 20383ab6a24SBarry Smith ierr = KSPSolve(red->ksp,red->xsub,red->ysub);CHKERRQ(ierr); 2043f457be1SHong Zhang ierr = VecResetArray(red->xsub);CHKERRQ(ierr); 2053f457be1SHong Zhang ierr = VecRestoreArray(red->xdup,&array);CHKERRQ(ierr); 2064b9ad928SBarry Smith 2073f457be1SHong Zhang /* place ysub's local array into ydup */ 2083f457be1SHong Zhang ierr = VecGetArray(red->ysub,&array);CHKERRQ(ierr); 2093f457be1SHong Zhang ierr = VecPlaceArray(red->ydup,(const PetscScalar*)array);CHKERRQ(ierr); 2103f457be1SHong Zhang 2113f457be1SHong Zhang /* scatter ydup to y */ 212ca9f406cSSatish Balay ierr = VecScatterBegin(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 213ca9f406cSSatish Balay ierr = VecScatterEnd(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2143f457be1SHong Zhang ierr = VecResetArray(red->ydup);CHKERRQ(ierr); 2153f457be1SHong Zhang ierr = VecRestoreArray(red->ysub,&array);CHKERRQ(ierr); 2164b9ad928SBarry Smith PetscFunctionReturn(0); 2174b9ad928SBarry Smith } 2184b9ad928SBarry Smith 219d88bfacbSStefano Zampini static PetscErrorCode PCApplyTranspose_Redundant(PC pc,Vec x,Vec y) 220d88bfacbSStefano Zampini { 221d88bfacbSStefano Zampini PC_Redundant *red = (PC_Redundant*)pc->data; 222d88bfacbSStefano Zampini PetscErrorCode ierr; 223d88bfacbSStefano Zampini PetscScalar *array; 224d88bfacbSStefano Zampini 225d88bfacbSStefano Zampini PetscFunctionBegin; 226d88bfacbSStefano Zampini if (!red->useparallelmat) { 227d88bfacbSStefano Zampini ierr = KSPSolveTranspose(red->ksp,x,y);CHKERRQ(ierr); 228d88bfacbSStefano Zampini PetscFunctionReturn(0); 229d88bfacbSStefano Zampini } 230d88bfacbSStefano Zampini 231d88bfacbSStefano Zampini /* scatter x to xdup */ 232d88bfacbSStefano Zampini ierr = VecScatterBegin(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 233d88bfacbSStefano Zampini ierr = VecScatterEnd(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 234d88bfacbSStefano Zampini 235d88bfacbSStefano Zampini /* place xdup's local array into xsub */ 236d88bfacbSStefano Zampini ierr = VecGetArray(red->xdup,&array);CHKERRQ(ierr); 237d88bfacbSStefano Zampini ierr = VecPlaceArray(red->xsub,(const PetscScalar*)array);CHKERRQ(ierr); 238d88bfacbSStefano Zampini 239d88bfacbSStefano Zampini /* apply preconditioner on each processor */ 240d88bfacbSStefano Zampini ierr = KSPSolveTranspose(red->ksp,red->xsub,red->ysub);CHKERRQ(ierr); 241d88bfacbSStefano Zampini ierr = VecResetArray(red->xsub);CHKERRQ(ierr); 242d88bfacbSStefano Zampini ierr = VecRestoreArray(red->xdup,&array);CHKERRQ(ierr); 243d88bfacbSStefano Zampini 244d88bfacbSStefano Zampini /* place ysub's local array into ydup */ 245d88bfacbSStefano Zampini ierr = VecGetArray(red->ysub,&array);CHKERRQ(ierr); 246d88bfacbSStefano Zampini ierr = VecPlaceArray(red->ydup,(const PetscScalar*)array);CHKERRQ(ierr); 247d88bfacbSStefano Zampini 248d88bfacbSStefano Zampini /* scatter ydup to y */ 249d88bfacbSStefano Zampini ierr = VecScatterBegin(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 250d88bfacbSStefano Zampini ierr = VecScatterEnd(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 251d88bfacbSStefano Zampini ierr = VecResetArray(red->ydup);CHKERRQ(ierr); 252d88bfacbSStefano Zampini ierr = VecRestoreArray(red->ysub,&array);CHKERRQ(ierr); 253d88bfacbSStefano Zampini PetscFunctionReturn(0); 254d88bfacbSStefano Zampini } 255d88bfacbSStefano Zampini 2561ea5a559SBarry Smith static PetscErrorCode PCReset_Redundant(PC pc) 2574b9ad928SBarry Smith { 2584b9ad928SBarry Smith PC_Redundant *red = (PC_Redundant*)pc->data; 259dfbe8321SBarry Smith PetscErrorCode ierr; 2604b9ad928SBarry Smith 2614b9ad928SBarry Smith PetscFunctionBegin; 2621b81debcSHong Zhang if (red->useparallelmat) { 2636bf464f9SBarry Smith ierr = VecScatterDestroy(&red->scatterin);CHKERRQ(ierr); 2646bf464f9SBarry Smith ierr = VecScatterDestroy(&red->scatterout);CHKERRQ(ierr); 2656bf464f9SBarry Smith ierr = VecDestroy(&red->ysub);CHKERRQ(ierr); 2666bf464f9SBarry Smith ierr = VecDestroy(&red->xsub);CHKERRQ(ierr); 2676bf464f9SBarry Smith ierr = VecDestroy(&red->xdup);CHKERRQ(ierr); 2686bf464f9SBarry Smith ierr = VecDestroy(&red->ydup);CHKERRQ(ierr); 2691b81debcSHong Zhang } 2706bf464f9SBarry Smith ierr = MatDestroy(&red->pmats);CHKERRQ(ierr); 2711b81debcSHong Zhang ierr = KSPReset(red->ksp);CHKERRQ(ierr); 2721ea5a559SBarry Smith PetscFunctionReturn(0); 2731ea5a559SBarry Smith } 2741ea5a559SBarry Smith 2751ea5a559SBarry Smith static PetscErrorCode PCDestroy_Redundant(PC pc) 2761ea5a559SBarry Smith { 2771ea5a559SBarry Smith PC_Redundant *red = (PC_Redundant*)pc->data; 2781ea5a559SBarry Smith PetscErrorCode ierr; 2791ea5a559SBarry Smith 2801ea5a559SBarry Smith PetscFunctionBegin; 2811ea5a559SBarry Smith ierr = PCReset_Redundant(pc);CHKERRQ(ierr); 2826bf464f9SBarry Smith ierr = KSPDestroy(&red->ksp);CHKERRQ(ierr); 2836bf464f9SBarry Smith ierr = PetscSubcommDestroy(&red->psubcomm);CHKERRQ(ierr); 284c31cb41cSBarry Smith ierr = PetscFree(pc->data);CHKERRQ(ierr); 2854b9ad928SBarry Smith PetscFunctionReturn(0); 2864b9ad928SBarry Smith } 2874b9ad928SBarry Smith 2884416b707SBarry Smith static PetscErrorCode PCSetFromOptions_Redundant(PetscOptionItems *PetscOptionsObject,PC pc) 2894b9ad928SBarry Smith { 290a98ce0f4SHong Zhang PetscErrorCode ierr; 291a98ce0f4SHong Zhang PC_Redundant *red = (PC_Redundant*)pc->data; 292a98ce0f4SHong Zhang 2934b9ad928SBarry Smith PetscFunctionBegin; 294e55864a3SBarry Smith ierr = PetscOptionsHead(PetscOptionsObject,"Redundant options");CHKERRQ(ierr); 29509a6bc64SHong Zhang ierr = PetscOptionsInt("-pc_redundant_number","Number of redundant pc","PCRedundantSetNumber",red->nsubcomm,&red->nsubcomm,0);CHKERRQ(ierr); 296a98ce0f4SHong Zhang ierr = PetscOptionsTail();CHKERRQ(ierr); 2974b9ad928SBarry Smith PetscFunctionReturn(0); 2984b9ad928SBarry Smith } 2994b9ad928SBarry Smith 300f7a08781SBarry Smith static PetscErrorCode PCRedundantSetNumber_Redundant(PC pc,PetscInt nreds) 30109a6bc64SHong Zhang { 30209a6bc64SHong Zhang PC_Redundant *red = (PC_Redundant*)pc->data; 30309a6bc64SHong Zhang 30409a6bc64SHong Zhang PetscFunctionBegin; 30509a6bc64SHong Zhang red->nsubcomm = nreds; 30609a6bc64SHong Zhang PetscFunctionReturn(0); 30709a6bc64SHong Zhang } 30809a6bc64SHong Zhang 30909a6bc64SHong Zhang /*@ 31009a6bc64SHong Zhang PCRedundantSetNumber - Sets the number of redundant preconditioner contexts. 31109a6bc64SHong Zhang 3123f9fe445SBarry Smith Logically Collective on PC 31309a6bc64SHong Zhang 31409a6bc64SHong Zhang Input Parameters: 31509a6bc64SHong Zhang + pc - the preconditioner context 3169b21d695SBarry Smith - nredundant - number of redundant preconditioner contexts; for example if you are using 64 MPI processes and 3179b21d695SBarry Smith use an nredundant of 4 there will be 4 parallel solves each on 16 = 64/4 processes. 31809a6bc64SHong Zhang 31909a6bc64SHong Zhang Level: advanced 32009a6bc64SHong Zhang 32109a6bc64SHong Zhang .keywords: PC, redundant solve 32209a6bc64SHong Zhang @*/ 3237087cfbeSBarry Smith PetscErrorCode PCRedundantSetNumber(PC pc,PetscInt nredundant) 32409a6bc64SHong Zhang { 3254ac538c5SBarry Smith PetscErrorCode ierr; 32609a6bc64SHong Zhang 32709a6bc64SHong Zhang PetscFunctionBegin; 3280700a824SBarry Smith PetscValidHeaderSpecific(pc,PC_CLASSID,1); 329ce94432eSBarry Smith if (nredundant <= 0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "num of redundant pc %D must be positive",nredundant); 3304ac538c5SBarry Smith ierr = PetscTryMethod(pc,"PCRedundantSetNumber_C",(PC,PetscInt),(pc,nredundant));CHKERRQ(ierr); 33109a6bc64SHong Zhang PetscFunctionReturn(0); 33209a6bc64SHong Zhang } 33309a6bc64SHong Zhang 334f7a08781SBarry Smith static PetscErrorCode PCRedundantSetScatter_Redundant(PC pc,VecScatter in,VecScatter out) 3354b9ad928SBarry Smith { 3364b9ad928SBarry Smith PC_Redundant *red = (PC_Redundant*)pc->data; 337dfbe8321SBarry Smith PetscErrorCode ierr; 3384b9ad928SBarry Smith 3394b9ad928SBarry Smith PetscFunctionBegin; 3404b9ad928SBarry Smith ierr = PetscObjectReference((PetscObject)in);CHKERRQ(ierr); 3416bf464f9SBarry Smith ierr = VecScatterDestroy(&red->scatterin);CHKERRQ(ierr); 3422fa5cd67SKarl Rupp 343c3122656SLisandro Dalcin red->scatterin = in; 3442fa5cd67SKarl Rupp 3454b9ad928SBarry Smith ierr = PetscObjectReference((PetscObject)out);CHKERRQ(ierr); 3466bf464f9SBarry Smith ierr = VecScatterDestroy(&red->scatterout);CHKERRQ(ierr); 347c3122656SLisandro Dalcin red->scatterout = out; 3484b9ad928SBarry Smith PetscFunctionReturn(0); 3494b9ad928SBarry Smith } 3504b9ad928SBarry Smith 3514b9ad928SBarry Smith /*@ 3524b9ad928SBarry Smith PCRedundantSetScatter - Sets the scatter used to copy values into the 3534b9ad928SBarry Smith redundant local solve and the scatter to move them back into the global 3544b9ad928SBarry Smith vector. 3554b9ad928SBarry Smith 3563f9fe445SBarry Smith Logically Collective on PC 3574b9ad928SBarry Smith 3584b9ad928SBarry Smith Input Parameters: 3594b9ad928SBarry Smith + pc - the preconditioner context 3604b9ad928SBarry Smith . in - the scatter to move the values in 3614b9ad928SBarry Smith - out - the scatter to move them out 3624b9ad928SBarry Smith 3634b9ad928SBarry Smith Level: advanced 3644b9ad928SBarry Smith 3654b9ad928SBarry Smith .keywords: PC, redundant solve 3664b9ad928SBarry Smith @*/ 3677087cfbeSBarry Smith PetscErrorCode PCRedundantSetScatter(PC pc,VecScatter in,VecScatter out) 3684b9ad928SBarry Smith { 3694ac538c5SBarry Smith PetscErrorCode ierr; 3704b9ad928SBarry Smith 3714b9ad928SBarry Smith PetscFunctionBegin; 3720700a824SBarry Smith PetscValidHeaderSpecific(pc,PC_CLASSID,1); 3730700a824SBarry Smith PetscValidHeaderSpecific(in,VEC_SCATTER_CLASSID,2); 3740700a824SBarry Smith PetscValidHeaderSpecific(out,VEC_SCATTER_CLASSID,3); 3754ac538c5SBarry Smith ierr = PetscTryMethod(pc,"PCRedundantSetScatter_C",(PC,VecScatter,VecScatter),(pc,in,out));CHKERRQ(ierr); 3764b9ad928SBarry Smith PetscFunctionReturn(0); 3774b9ad928SBarry Smith } 3784b9ad928SBarry Smith 379f7a08781SBarry Smith static PetscErrorCode PCRedundantGetKSP_Redundant(PC pc,KSP *innerksp) 3804b9ad928SBarry Smith { 3815f06b7aaSBarry Smith PetscErrorCode ierr; 3824b9ad928SBarry Smith PC_Redundant *red = (PC_Redundant*)pc->data; 38375024027SHong Zhang MPI_Comm comm,subcomm; 38475024027SHong Zhang const char *prefix; 3854b9ad928SBarry Smith 3864b9ad928SBarry Smith PetscFunctionBegin; 38775024027SHong Zhang if (!red->psubcomm) { 388e5acf8a4SHong Zhang ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr); 389e5acf8a4SHong Zhang 39075024027SHong Zhang ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 39175024027SHong Zhang ierr = PetscSubcommCreate(comm,&red->psubcomm);CHKERRQ(ierr); 39275024027SHong Zhang ierr = PetscSubcommSetNumber(red->psubcomm,red->nsubcomm);CHKERRQ(ierr); 39375024027SHong Zhang ierr = PetscSubcommSetType(red->psubcomm,PETSC_SUBCOMM_CONTIGUOUS);CHKERRQ(ierr); 394e5acf8a4SHong Zhang 395e5acf8a4SHong Zhang ierr = PetscSubcommSetOptionsPrefix(red->psubcomm,prefix);CHKERRQ(ierr); 396e5acf8a4SHong Zhang ierr = PetscSubcommSetFromOptions(red->psubcomm);CHKERRQ(ierr); 39775024027SHong Zhang ierr = PetscLogObjectMemory((PetscObject)pc,sizeof(PetscSubcomm));CHKERRQ(ierr); 39875024027SHong Zhang 39975024027SHong Zhang /* create a new PC that processors in each subcomm have copy of */ 40075024027SHong Zhang subcomm = PetscSubcommChild(red->psubcomm); 40175024027SHong Zhang 40275024027SHong Zhang ierr = KSPCreate(subcomm,&red->ksp);CHKERRQ(ierr); 40375024027SHong Zhang ierr = KSPSetErrorIfNotConverged(red->ksp,pc->erroriffailure);CHKERRQ(ierr); 40475024027SHong Zhang ierr = PetscObjectIncrementTabLevel((PetscObject)red->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 40575024027SHong Zhang ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)red->ksp);CHKERRQ(ierr); 40675024027SHong Zhang ierr = KSPSetType(red->ksp,KSPPREONLY);CHKERRQ(ierr); 40775024027SHong Zhang ierr = KSPGetPC(red->ksp,&red->pc);CHKERRQ(ierr); 40875024027SHong Zhang ierr = PCSetType(red->pc,PCLU);CHKERRQ(ierr); 409753b7fb9SBarry Smith if (red->shifttypeset) { 410753b7fb9SBarry Smith ierr = PCFactorSetShiftType(red->pc,red->shifttype);CHKERRQ(ierr); 411753b7fb9SBarry Smith red->shifttypeset = PETSC_FALSE; 412753b7fb9SBarry Smith } 41375024027SHong Zhang ierr = KSPSetOptionsPrefix(red->ksp,prefix);CHKERRQ(ierr); 41475024027SHong Zhang ierr = KSPAppendOptionsPrefix(red->ksp,"redundant_");CHKERRQ(ierr); 41575024027SHong Zhang } 41683ab6a24SBarry Smith *innerksp = red->ksp; 4174b9ad928SBarry Smith PetscFunctionReturn(0); 4184b9ad928SBarry Smith } 4194b9ad928SBarry Smith 4204b9ad928SBarry Smith /*@ 42183ab6a24SBarry Smith PCRedundantGetKSP - Gets the less parallel KSP created by the redundant PC. 4224b9ad928SBarry Smith 4234b9ad928SBarry Smith Not Collective 4244b9ad928SBarry Smith 4254b9ad928SBarry Smith Input Parameter: 4264b9ad928SBarry Smith . pc - the preconditioner context 4274b9ad928SBarry Smith 4284b9ad928SBarry Smith Output Parameter: 42983ab6a24SBarry Smith . innerksp - the KSP on the smaller set of processes 4304b9ad928SBarry Smith 4314b9ad928SBarry Smith Level: advanced 4324b9ad928SBarry Smith 4334b9ad928SBarry Smith .keywords: PC, redundant solve 4344b9ad928SBarry Smith @*/ 43583ab6a24SBarry Smith PetscErrorCode PCRedundantGetKSP(PC pc,KSP *innerksp) 4364b9ad928SBarry Smith { 4374ac538c5SBarry Smith PetscErrorCode ierr; 4384b9ad928SBarry Smith 4394b9ad928SBarry Smith PetscFunctionBegin; 4400700a824SBarry Smith PetscValidHeaderSpecific(pc,PC_CLASSID,1); 44183ab6a24SBarry Smith PetscValidPointer(innerksp,2); 442753b7fb9SBarry Smith ierr = PetscUseMethod(pc,"PCRedundantGetKSP_C",(PC,KSP*),(pc,innerksp));CHKERRQ(ierr); 4434b9ad928SBarry Smith PetscFunctionReturn(0); 4444b9ad928SBarry Smith } 4454b9ad928SBarry Smith 446f7a08781SBarry Smith static PetscErrorCode PCRedundantGetOperators_Redundant(PC pc,Mat *mat,Mat *pmat) 4474b9ad928SBarry Smith { 4484b9ad928SBarry Smith PC_Redundant *red = (PC_Redundant*)pc->data; 4494b9ad928SBarry Smith 4504b9ad928SBarry Smith PetscFunctionBegin; 451b3804887SHong Zhang if (mat) *mat = red->pmats; 452b3804887SHong Zhang if (pmat) *pmat = red->pmats; 4534b9ad928SBarry Smith PetscFunctionReturn(0); 4544b9ad928SBarry Smith } 4554b9ad928SBarry Smith 4564b9ad928SBarry Smith /*@ 4574b9ad928SBarry Smith PCRedundantGetOperators - gets the sequential matrix and preconditioner matrix 4584b9ad928SBarry Smith 4594b9ad928SBarry Smith Not Collective 4604b9ad928SBarry Smith 4614b9ad928SBarry Smith Input Parameter: 4624b9ad928SBarry Smith . pc - the preconditioner context 4634b9ad928SBarry Smith 4644b9ad928SBarry Smith Output Parameters: 4654b9ad928SBarry Smith + mat - the matrix 4664b9ad928SBarry Smith - pmat - the (possibly different) preconditioner matrix 4674b9ad928SBarry Smith 4684b9ad928SBarry Smith Level: advanced 4694b9ad928SBarry Smith 4704b9ad928SBarry Smith .keywords: PC, redundant solve 4714b9ad928SBarry Smith @*/ 4727087cfbeSBarry Smith PetscErrorCode PCRedundantGetOperators(PC pc,Mat *mat,Mat *pmat) 4734b9ad928SBarry Smith { 4744ac538c5SBarry Smith PetscErrorCode ierr; 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); 480163d334eSBarry Smith ierr = PetscUseMethod(pc,"PCRedundantGetOperators_C",(PC,Mat*,Mat*),(pc,mat,pmat));CHKERRQ(ierr); 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 49683ab6a24SBarry Smith Notes: The default KSP is preonly and the default PC is LU. 49783ab6a24SBarry Smith 498753b7fb9SBarry Smith PCFactorSetShiftType() applied to this PC will convey they shift type into the inner PC if it is factorization based. 499753b7fb9SBarry Smith 50083ab6a24SBarry Smith Developer Notes: Note that PCSetInitialGuessNonzero() is not used by this class but likely should be. 5019cfaa89bSBarry Smith 50237a17b4dSBarry Smith .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCRedundantSetScatter(), 50383ab6a24SBarry Smith PCRedundantGetKSP(), PCRedundantGetOperators(), PCRedundantSetNumber() 50437a17b4dSBarry Smith M*/ 50537a17b4dSBarry Smith 5068cc058d9SJed Brown PETSC_EXTERN PetscErrorCode PCCreate_Redundant(PC pc) 5074b9ad928SBarry Smith { 508dfbe8321SBarry Smith PetscErrorCode ierr; 5094b9ad928SBarry Smith PC_Redundant *red; 51069db28dcSHong Zhang PetscMPIInt size; 5113f457be1SHong Zhang 5124b9ad928SBarry Smith PetscFunctionBegin; 513b00a9115SJed Brown ierr = PetscNewLog(pc,&red);CHKERRQ(ierr); 514ce94432eSBarry Smith ierr = MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);CHKERRQ(ierr); 5152fa5cd67SKarl Rupp 51669db28dcSHong Zhang red->nsubcomm = size; 5174b9ad928SBarry Smith red->useparallelmat = PETSC_TRUE; 5181fbd8f88SHong Zhang pc->data = (void*)red; 5194b9ad928SBarry Smith 5204b9ad928SBarry Smith pc->ops->apply = PCApply_Redundant; 521d88bfacbSStefano Zampini pc->ops->applytranspose = PCApplyTranspose_Redundant; 5224b9ad928SBarry Smith pc->ops->setup = PCSetUp_Redundant; 5234b9ad928SBarry Smith pc->ops->destroy = PCDestroy_Redundant; 5241ea5a559SBarry Smith pc->ops->reset = PCReset_Redundant; 5254b9ad928SBarry Smith pc->ops->setfromoptions = PCSetFromOptions_Redundant; 5264b9ad928SBarry Smith pc->ops->view = PCView_Redundant; 5272fa5cd67SKarl Rupp 528bdf89e91SBarry Smith ierr = PetscObjectComposeFunction((PetscObject)pc,"PCRedundantSetScatter_C",PCRedundantSetScatter_Redundant);CHKERRQ(ierr); 529bdf89e91SBarry Smith ierr = PetscObjectComposeFunction((PetscObject)pc,"PCRedundantSetNumber_C",PCRedundantSetNumber_Redundant);CHKERRQ(ierr); 530bdf89e91SBarry Smith ierr = PetscObjectComposeFunction((PetscObject)pc,"PCRedundantGetKSP_C",PCRedundantGetKSP_Redundant);CHKERRQ(ierr); 531bdf89e91SBarry Smith ierr = PetscObjectComposeFunction((PetscObject)pc,"PCRedundantGetOperators_C",PCRedundantGetOperators_Redundant);CHKERRQ(ierr); 532753b7fb9SBarry Smith ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftType_C",PCFactorSetShiftType_Redundant);CHKERRQ(ierr); 5334b9ad928SBarry Smith PetscFunctionReturn(0); 5344b9ad928SBarry Smith } 535b2573a8aSBarry Smith 536