xref: /petsc/src/ksp/pc/impls/redundant/redundant.c (revision c0decd05c6848b80907752eef350b55c8c90e696)
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) {
51efd4aadfSBarry Smith       ierr = PetscViewerASCIIPrintf(viewer,"  Not yet setup\n");CHKERRQ(ierr);
5203ccd0b4SBarry Smith     } else {
53efd4aadfSBarry 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 
6819b3b6edSHong Zhang #include <../src/mat/impls/aij/mpi/mpiaij.h>
696849ba73SBarry Smith static PetscErrorCode PCSetUp_Redundant(PC pc)
704b9ad928SBarry Smith {
714b9ad928SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
72dfbe8321SBarry Smith   PetscErrorCode ierr;
731b81debcSHong Zhang   PetscInt       mstart,mend,mlocal,M;
7413f74950SBarry Smith   PetscMPIInt    size;
75ce94432eSBarry Smith   MPI_Comm       comm,subcomm;
76ddc54837SHong Zhang   Vec            x;
773f457be1SHong Zhang 
784b9ad928SBarry Smith   PetscFunctionBegin;
79ce94432eSBarry Smith   ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
80ddc54837SHong Zhang 
81ddc54837SHong Zhang   /* if pmatrix set by user is sequential then we do not need to gather the parallel matrix */
82ddc54837SHong Zhang   ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
83ddc54837SHong Zhang   if (size == 1) red->useparallelmat = PETSC_FALSE;
841fbd8f88SHong Zhang 
854b9ad928SBarry Smith   if (!pc->setupcalled) {
861b81debcSHong Zhang     PetscInt mloc_sub;
8775024027SHong Zhang     if (!red->psubcomm) { /* create red->psubcomm, new ksp and pc over subcomm */
8875024027SHong Zhang       KSP ksp;
8975024027SHong Zhang       ierr = PCRedundantGetKSP(pc,&ksp);CHKERRQ(ierr);
901b81debcSHong Zhang     }
9175024027SHong Zhang     subcomm = PetscSubcommChild(red->psubcomm);
921fbd8f88SHong Zhang 
931b81debcSHong Zhang     if (red->useparallelmat) {
941b81debcSHong Zhang       /* grab the parallel matrix and put it into processors of a subcomminicator */
9553cd1579SHong Zhang       ierr = MatCreateRedundantMatrix(pc->pmat,red->psubcomm->n,subcomm,MAT_INITIAL_MATRIX,&red->pmats);CHKERRQ(ierr);
96b85f2e9bSHong Zhang 
97b85f2e9bSHong Zhang       ierr = MPI_Comm_size(subcomm,&size);CHKERRQ(ierr);
98b85f2e9bSHong Zhang       if (size > 1) {
99b85f2e9bSHong Zhang         PetscBool foundpack;
100b85f2e9bSHong Zhang         ierr = MatGetFactorAvailable(red->pmats,NULL,MAT_FACTOR_LU,&foundpack);CHKERRQ(ierr);
101b85f2e9bSHong Zhang         if (!foundpack) { /* reset default ksp and pc */
102b85f2e9bSHong Zhang           ierr = KSPSetType(red->ksp,KSPGMRES);CHKERRQ(ierr);
103b85f2e9bSHong Zhang           ierr = PCSetType(red->pc,PCBJACOBI);CHKERRQ(ierr);
104c1619fb6SBarry Smith         } else {
1053ca39a21SBarry Smith           ierr = PCFactorSetMatSolverType(red->pc,NULL);CHKERRQ(ierr);
106b85f2e9bSHong Zhang         }
107b85f2e9bSHong Zhang       }
108b85f2e9bSHong Zhang 
10923ee1639SBarry Smith       ierr = KSPSetOperators(red->ksp,red->pmats,red->pmats);CHKERRQ(ierr);
1104b9ad928SBarry Smith 
1111b81debcSHong Zhang       /* get working vectors xsub and ysub */
1122a7a6963SBarry Smith       ierr = MatCreateVecs(red->pmats,&red->xsub,&red->ysub);CHKERRQ(ierr);
1132fa5cd67SKarl Rupp 
1148b96b0d2SHong Zhang       /* create working vectors xdup and ydup.
1158b96b0d2SHong Zhang        xdup concatenates all xsub's contigously to form a mpi vector over dupcomm  (see PetscSubcommCreate_interlaced())
1168b96b0d2SHong Zhang        ydup concatenates all ysub and has empty local arrays because ysub's arrays will be place into it.
117ce94432eSBarry Smith        Note: we use communicator dupcomm, not PetscObjectComm((PetscObject)pc)! */
1181b81debcSHong Zhang       ierr = MatGetLocalSize(red->pmats,&mloc_sub,NULL);CHKERRQ(ierr);
11936be1a5eSBarry Smith       ierr = VecCreateMPI(PetscSubcommContiguousParent(red->psubcomm),mloc_sub,PETSC_DECIDE,&red->xdup);CHKERRQ(ierr);
12036be1a5eSBarry Smith       ierr = VecCreateMPIWithArray(PetscSubcommContiguousParent(red->psubcomm),1,mloc_sub,PETSC_DECIDE,NULL,&red->ydup);CHKERRQ(ierr);
1213f457be1SHong Zhang 
122f68be91cSHong Zhang       /* create vecscatters */
123f68be91cSHong Zhang       if (!red->scatterin) { /* efficiency of scatterin is independent from psubcomm_type! */
1243f457be1SHong Zhang         IS       is1,is2;
1253f457be1SHong Zhang         PetscInt *idx1,*idx2,i,j,k;
12645fc02eaSBarry Smith 
1272a7a6963SBarry Smith         ierr = MatCreateVecs(pc->pmat,&x,0);CHKERRQ(ierr);
1281b81debcSHong Zhang         ierr = VecGetSize(x,&M);CHKERRQ(ierr);
1291b81debcSHong Zhang         ierr = VecGetOwnershipRange(x,&mstart,&mend);CHKERRQ(ierr);
1301b81debcSHong Zhang         mlocal = mend - mstart;
131dcca6d9dSJed Brown         ierr = PetscMalloc2(red->psubcomm->n*mlocal,&idx1,red->psubcomm->n*mlocal,&idx2);CHKERRQ(ierr);
1323f457be1SHong Zhang         j    = 0;
1331fbd8f88SHong Zhang         for (k=0; k<red->psubcomm->n; k++) {
1343f457be1SHong Zhang           for (i=mstart; i<mend; i++) {
1353f457be1SHong Zhang             idx1[j]   = i;
136ddc54837SHong Zhang             idx2[j++] = i + M*k;
1373f457be1SHong Zhang           }
1383f457be1SHong Zhang         }
13970b3c8c7SBarry Smith         ierr = ISCreateGeneral(comm,red->psubcomm->n*mlocal,idx1,PETSC_COPY_VALUES,&is1);CHKERRQ(ierr);
14070b3c8c7SBarry Smith         ierr = ISCreateGeneral(comm,red->psubcomm->n*mlocal,idx2,PETSC_COPY_VALUES,&is2);CHKERRQ(ierr);
14135928de7SBarry Smith         ierr = VecScatterCreateWithData(x,is1,red->xdup,is2,&red->scatterin);CHKERRQ(ierr);
142fcfd50ebSBarry Smith         ierr = ISDestroy(&is1);CHKERRQ(ierr);
143fcfd50ebSBarry Smith         ierr = ISDestroy(&is2);CHKERRQ(ierr);
1443f457be1SHong Zhang 
1456909748bSHong Zhang         /* Impl below is good for PETSC_SUBCOMM_INTERLACED (no inter-process communication) and PETSC_SUBCOMM_CONTIGUOUS (communication within subcomm) */
146ddc54837SHong Zhang         ierr = ISCreateStride(comm,mlocal,mstart+ red->psubcomm->color*M,1,&is1);CHKERRQ(ierr);
1473f457be1SHong Zhang         ierr = ISCreateStride(comm,mlocal,mstart,1,&is2);CHKERRQ(ierr);
14835928de7SBarry Smith         ierr = VecScatterCreateWithData(red->xdup,is1,x,is2,&red->scatterout);CHKERRQ(ierr);
149fcfd50ebSBarry Smith         ierr = ISDestroy(&is1);CHKERRQ(ierr);
150fcfd50ebSBarry Smith         ierr = ISDestroy(&is2);CHKERRQ(ierr);
1511d79065fSBarry Smith         ierr = PetscFree2(idx1,idx2);CHKERRQ(ierr);
152ddc54837SHong Zhang         ierr = VecDestroy(&x);CHKERRQ(ierr);
1531b81debcSHong Zhang       }
154ab661555SHong Zhang     } else { /* !red->useparallelmat */
15523ee1639SBarry Smith       ierr = KSPSetOperators(red->ksp,pc->mat,pc->pmat);CHKERRQ(ierr);
1561b81debcSHong Zhang     }
157ab661555SHong Zhang   } else { /* pc->setupcalled */
1584b9ad928SBarry Smith     if (red->useparallelmat) {
159ab661555SHong Zhang       MatReuse       reuse;
1601b81debcSHong Zhang       /* grab the parallel matrix and put it into processors of a subcomminicator */
1611b81debcSHong Zhang       /*--------------------------------------------------------------------------*/
162ab661555SHong Zhang       if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1634b9ad928SBarry Smith         /* destroy old matrices */
1646bf464f9SBarry Smith         ierr  = MatDestroy(&red->pmats);CHKERRQ(ierr);
165ab661555SHong Zhang         reuse = MAT_INITIAL_MATRIX;
1664b9ad928SBarry Smith       } else {
167ab661555SHong Zhang         reuse = MAT_REUSE_MATRIX;
168ab661555SHong Zhang       }
169306c2d5bSBarry Smith       ierr = MatCreateRedundantMatrix(pc->pmat,red->psubcomm->n,PetscSubcommChild(red->psubcomm),reuse,&red->pmats);CHKERRQ(ierr);
17023ee1639SBarry Smith       ierr = KSPSetOperators(red->ksp,red->pmats,red->pmats);CHKERRQ(ierr);
171ab661555SHong Zhang     } else { /* !red->useparallelmat */
17223ee1639SBarry Smith       ierr = KSPSetOperators(red->ksp,pc->mat,pc->pmat);CHKERRQ(ierr);
1734b9ad928SBarry Smith     }
174ab661555SHong Zhang   }
1751b81debcSHong Zhang 
1760c24e6a1SHong Zhang   if (pc->setfromoptionscalled) {
1773e065800SHong Zhang     ierr = KSPSetFromOptions(red->ksp);CHKERRQ(ierr);
1780c24e6a1SHong Zhang   }
1793e065800SHong Zhang   ierr = KSPSetUp(red->ksp);CHKERRQ(ierr);
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;
186dfbe8321SBarry Smith   PetscErrorCode ierr;
1873f457be1SHong Zhang   PetscScalar    *array;
1884b9ad928SBarry Smith 
1894b9ad928SBarry Smith   PetscFunctionBegin;
190ddc54837SHong Zhang   if (!red->useparallelmat) {
191ddc54837SHong Zhang     ierr = KSPSolve(red->ksp,x,y);CHKERRQ(ierr);
192*c0decd05SBarry Smith     ierr = KSPCheckSolve(red->ksp,pc,y);CHKERRQ(ierr);
193ddc54837SHong Zhang     PetscFunctionReturn(0);
194ddc54837SHong Zhang   }
195ddc54837SHong Zhang 
1963f457be1SHong Zhang   /* scatter x to xdup */
197ca9f406cSSatish Balay   ierr = VecScatterBegin(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
198ca9f406cSSatish Balay   ierr = VecScatterEnd(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1993f457be1SHong Zhang 
2003f457be1SHong Zhang   /* place xdup's local array into xsub */
2013f457be1SHong Zhang   ierr = VecGetArray(red->xdup,&array);CHKERRQ(ierr);
2023f457be1SHong Zhang   ierr = VecPlaceArray(red->xsub,(const PetscScalar*)array);CHKERRQ(ierr);
2034b9ad928SBarry Smith 
2044b9ad928SBarry Smith   /* apply preconditioner on each processor */
20583ab6a24SBarry Smith   ierr = KSPSolve(red->ksp,red->xsub,red->ysub);CHKERRQ(ierr);
206*c0decd05SBarry Smith   ierr = KSPCheckSolve(red->ksp,pc,red->ysub);CHKERRQ(ierr);
2073f457be1SHong Zhang   ierr = VecResetArray(red->xsub);CHKERRQ(ierr);
2083f457be1SHong Zhang   ierr = VecRestoreArray(red->xdup,&array);CHKERRQ(ierr);
2094b9ad928SBarry Smith 
2103f457be1SHong Zhang   /* place ysub's local array into ydup */
2113f457be1SHong Zhang   ierr = VecGetArray(red->ysub,&array);CHKERRQ(ierr);
2123f457be1SHong Zhang   ierr = VecPlaceArray(red->ydup,(const PetscScalar*)array);CHKERRQ(ierr);
2133f457be1SHong Zhang 
2143f457be1SHong Zhang   /* scatter ydup to y */
215ca9f406cSSatish Balay   ierr = VecScatterBegin(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
216ca9f406cSSatish Balay   ierr = VecScatterEnd(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
2173f457be1SHong Zhang   ierr = VecResetArray(red->ydup);CHKERRQ(ierr);
2183f457be1SHong Zhang   ierr = VecRestoreArray(red->ysub,&array);CHKERRQ(ierr);
2194b9ad928SBarry Smith   PetscFunctionReturn(0);
2204b9ad928SBarry Smith }
2214b9ad928SBarry Smith 
222d88bfacbSStefano Zampini static PetscErrorCode PCApplyTranspose_Redundant(PC pc,Vec x,Vec y)
223d88bfacbSStefano Zampini {
224d88bfacbSStefano Zampini   PC_Redundant   *red = (PC_Redundant*)pc->data;
225d88bfacbSStefano Zampini   PetscErrorCode ierr;
226d88bfacbSStefano Zampini   PetscScalar    *array;
227d88bfacbSStefano Zampini 
228d88bfacbSStefano Zampini   PetscFunctionBegin;
229d88bfacbSStefano Zampini   if (!red->useparallelmat) {
230d88bfacbSStefano Zampini     ierr = KSPSolveTranspose(red->ksp,x,y);CHKERRQ(ierr);
231*c0decd05SBarry Smith     ierr = KSPCheckSolve(red->ksp,pc,y);CHKERRQ(ierr);
232d88bfacbSStefano Zampini     PetscFunctionReturn(0);
233d88bfacbSStefano Zampini   }
234d88bfacbSStefano Zampini 
235d88bfacbSStefano Zampini   /* scatter x to xdup */
236d88bfacbSStefano Zampini   ierr = VecScatterBegin(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
237d88bfacbSStefano Zampini   ierr = VecScatterEnd(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
238d88bfacbSStefano Zampini 
239d88bfacbSStefano Zampini   /* place xdup's local array into xsub */
240d88bfacbSStefano Zampini   ierr = VecGetArray(red->xdup,&array);CHKERRQ(ierr);
241d88bfacbSStefano Zampini   ierr = VecPlaceArray(red->xsub,(const PetscScalar*)array);CHKERRQ(ierr);
242d88bfacbSStefano Zampini 
243d88bfacbSStefano Zampini   /* apply preconditioner on each processor */
244d88bfacbSStefano Zampini   ierr = KSPSolveTranspose(red->ksp,red->xsub,red->ysub);CHKERRQ(ierr);
245*c0decd05SBarry Smith   ierr = KSPCheckSolve(red->ksp,pc,red->ysub);CHKERRQ(ierr);
246d88bfacbSStefano Zampini   ierr = VecResetArray(red->xsub);CHKERRQ(ierr);
247d88bfacbSStefano Zampini   ierr = VecRestoreArray(red->xdup,&array);CHKERRQ(ierr);
248d88bfacbSStefano Zampini 
249d88bfacbSStefano Zampini   /* place ysub's local array into ydup */
250d88bfacbSStefano Zampini   ierr = VecGetArray(red->ysub,&array);CHKERRQ(ierr);
251d88bfacbSStefano Zampini   ierr = VecPlaceArray(red->ydup,(const PetscScalar*)array);CHKERRQ(ierr);
252d88bfacbSStefano Zampini 
253d88bfacbSStefano Zampini   /* scatter ydup to y */
254d88bfacbSStefano Zampini   ierr = VecScatterBegin(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
255d88bfacbSStefano Zampini   ierr = VecScatterEnd(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
256d88bfacbSStefano Zampini   ierr = VecResetArray(red->ydup);CHKERRQ(ierr);
257d88bfacbSStefano Zampini   ierr = VecRestoreArray(red->ysub,&array);CHKERRQ(ierr);
258d88bfacbSStefano Zampini   PetscFunctionReturn(0);
259d88bfacbSStefano Zampini }
260d88bfacbSStefano Zampini 
2611ea5a559SBarry Smith static PetscErrorCode PCReset_Redundant(PC pc)
2624b9ad928SBarry Smith {
2634b9ad928SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
264dfbe8321SBarry Smith   PetscErrorCode ierr;
2654b9ad928SBarry Smith 
2664b9ad928SBarry Smith   PetscFunctionBegin;
2671b81debcSHong Zhang   if (red->useparallelmat) {
2686bf464f9SBarry Smith     ierr = VecScatterDestroy(&red->scatterin);CHKERRQ(ierr);
2696bf464f9SBarry Smith     ierr = VecScatterDestroy(&red->scatterout);CHKERRQ(ierr);
2706bf464f9SBarry Smith     ierr = VecDestroy(&red->ysub);CHKERRQ(ierr);
2716bf464f9SBarry Smith     ierr = VecDestroy(&red->xsub);CHKERRQ(ierr);
2726bf464f9SBarry Smith     ierr = VecDestroy(&red->xdup);CHKERRQ(ierr);
2736bf464f9SBarry Smith     ierr = VecDestroy(&red->ydup);CHKERRQ(ierr);
2741b81debcSHong Zhang   }
2756bf464f9SBarry Smith   ierr = MatDestroy(&red->pmats);CHKERRQ(ierr);
2761b81debcSHong Zhang   ierr = KSPReset(red->ksp);CHKERRQ(ierr);
2771ea5a559SBarry Smith   PetscFunctionReturn(0);
2781ea5a559SBarry Smith }
2791ea5a559SBarry Smith 
2801ea5a559SBarry Smith static PetscErrorCode PCDestroy_Redundant(PC pc)
2811ea5a559SBarry Smith {
2821ea5a559SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
2831ea5a559SBarry Smith   PetscErrorCode ierr;
2841ea5a559SBarry Smith 
2851ea5a559SBarry Smith   PetscFunctionBegin;
2861ea5a559SBarry Smith   ierr = PCReset_Redundant(pc);CHKERRQ(ierr);
2876bf464f9SBarry Smith   ierr = KSPDestroy(&red->ksp);CHKERRQ(ierr);
2886bf464f9SBarry Smith   ierr = PetscSubcommDestroy(&red->psubcomm);CHKERRQ(ierr);
289c31cb41cSBarry Smith   ierr = PetscFree(pc->data);CHKERRQ(ierr);
2904b9ad928SBarry Smith   PetscFunctionReturn(0);
2914b9ad928SBarry Smith }
2924b9ad928SBarry Smith 
2934416b707SBarry Smith static PetscErrorCode PCSetFromOptions_Redundant(PetscOptionItems *PetscOptionsObject,PC pc)
2944b9ad928SBarry Smith {
295a98ce0f4SHong Zhang   PetscErrorCode ierr;
296a98ce0f4SHong Zhang   PC_Redundant   *red = (PC_Redundant*)pc->data;
297a98ce0f4SHong Zhang 
2984b9ad928SBarry Smith   PetscFunctionBegin;
299e55864a3SBarry Smith   ierr = PetscOptionsHead(PetscOptionsObject,"Redundant options");CHKERRQ(ierr);
30009a6bc64SHong Zhang   ierr = PetscOptionsInt("-pc_redundant_number","Number of redundant pc","PCRedundantSetNumber",red->nsubcomm,&red->nsubcomm,0);CHKERRQ(ierr);
301a98ce0f4SHong Zhang   ierr = PetscOptionsTail();CHKERRQ(ierr);
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 .keywords: PC, redundant solve
32709a6bc64SHong Zhang @*/
3287087cfbeSBarry Smith PetscErrorCode PCRedundantSetNumber(PC pc,PetscInt nredundant)
32909a6bc64SHong Zhang {
3304ac538c5SBarry Smith   PetscErrorCode ierr;
33109a6bc64SHong Zhang 
33209a6bc64SHong Zhang   PetscFunctionBegin;
3330700a824SBarry Smith   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
334ce94432eSBarry Smith   if (nredundant <= 0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "num of redundant pc %D must be positive",nredundant);
3354ac538c5SBarry Smith   ierr = PetscTryMethod(pc,"PCRedundantSetNumber_C",(PC,PetscInt),(pc,nredundant));CHKERRQ(ierr);
33609a6bc64SHong Zhang   PetscFunctionReturn(0);
33709a6bc64SHong Zhang }
33809a6bc64SHong Zhang 
339f7a08781SBarry Smith static PetscErrorCode PCRedundantSetScatter_Redundant(PC pc,VecScatter in,VecScatter out)
3404b9ad928SBarry Smith {
3414b9ad928SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
342dfbe8321SBarry Smith   PetscErrorCode ierr;
3434b9ad928SBarry Smith 
3444b9ad928SBarry Smith   PetscFunctionBegin;
3454b9ad928SBarry Smith   ierr = PetscObjectReference((PetscObject)in);CHKERRQ(ierr);
3466bf464f9SBarry Smith   ierr = VecScatterDestroy(&red->scatterin);CHKERRQ(ierr);
3472fa5cd67SKarl Rupp 
348c3122656SLisandro Dalcin   red->scatterin  = in;
3492fa5cd67SKarl Rupp 
3504b9ad928SBarry Smith   ierr = PetscObjectReference((PetscObject)out);CHKERRQ(ierr);
3516bf464f9SBarry Smith   ierr = VecScatterDestroy(&red->scatterout);CHKERRQ(ierr);
352c3122656SLisandro Dalcin   red->scatterout = out;
3534b9ad928SBarry Smith   PetscFunctionReturn(0);
3544b9ad928SBarry Smith }
3554b9ad928SBarry Smith 
3564b9ad928SBarry Smith /*@
3574b9ad928SBarry Smith    PCRedundantSetScatter - Sets the scatter used to copy values into the
3584b9ad928SBarry Smith      redundant local solve and the scatter to move them back into the global
3594b9ad928SBarry Smith      vector.
3604b9ad928SBarry Smith 
3613f9fe445SBarry Smith    Logically Collective on PC
3624b9ad928SBarry Smith 
3634b9ad928SBarry Smith    Input Parameters:
3644b9ad928SBarry Smith +  pc - the preconditioner context
3654b9ad928SBarry Smith .  in - the scatter to move the values in
3664b9ad928SBarry Smith -  out - the scatter to move them out
3674b9ad928SBarry Smith 
3684b9ad928SBarry Smith    Level: advanced
3694b9ad928SBarry Smith 
3704b9ad928SBarry Smith .keywords: PC, redundant solve
3714b9ad928SBarry Smith @*/
3727087cfbeSBarry Smith PetscErrorCode PCRedundantSetScatter(PC pc,VecScatter in,VecScatter out)
3734b9ad928SBarry Smith {
3744ac538c5SBarry Smith   PetscErrorCode ierr;
3754b9ad928SBarry Smith 
3764b9ad928SBarry Smith   PetscFunctionBegin;
3770700a824SBarry Smith   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
3780700a824SBarry Smith   PetscValidHeaderSpecific(in,VEC_SCATTER_CLASSID,2);
3790700a824SBarry Smith   PetscValidHeaderSpecific(out,VEC_SCATTER_CLASSID,3);
3804ac538c5SBarry Smith   ierr = PetscTryMethod(pc,"PCRedundantSetScatter_C",(PC,VecScatter,VecScatter),(pc,in,out));CHKERRQ(ierr);
3814b9ad928SBarry Smith   PetscFunctionReturn(0);
3824b9ad928SBarry Smith }
3834b9ad928SBarry Smith 
384f7a08781SBarry Smith static PetscErrorCode PCRedundantGetKSP_Redundant(PC pc,KSP *innerksp)
3854b9ad928SBarry Smith {
3865f06b7aaSBarry Smith   PetscErrorCode ierr;
3874b9ad928SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
38875024027SHong Zhang   MPI_Comm       comm,subcomm;
38975024027SHong Zhang   const char     *prefix;
3904b9ad928SBarry Smith 
3914b9ad928SBarry Smith   PetscFunctionBegin;
39275024027SHong Zhang   if (!red->psubcomm) {
393e5acf8a4SHong Zhang     ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr);
394e5acf8a4SHong Zhang 
39575024027SHong Zhang     ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
39675024027SHong Zhang     ierr = PetscSubcommCreate(comm,&red->psubcomm);CHKERRQ(ierr);
39775024027SHong Zhang     ierr = PetscSubcommSetNumber(red->psubcomm,red->nsubcomm);CHKERRQ(ierr);
39875024027SHong Zhang     ierr = PetscSubcommSetType(red->psubcomm,PETSC_SUBCOMM_CONTIGUOUS);CHKERRQ(ierr);
399e5acf8a4SHong Zhang 
400e5acf8a4SHong Zhang     ierr = PetscSubcommSetOptionsPrefix(red->psubcomm,prefix);CHKERRQ(ierr);
401e5acf8a4SHong Zhang     ierr = PetscSubcommSetFromOptions(red->psubcomm);CHKERRQ(ierr);
40275024027SHong Zhang     ierr = PetscLogObjectMemory((PetscObject)pc,sizeof(PetscSubcomm));CHKERRQ(ierr);
40375024027SHong Zhang 
40475024027SHong Zhang     /* create a new PC that processors in each subcomm have copy of */
40575024027SHong Zhang     subcomm = PetscSubcommChild(red->psubcomm);
40675024027SHong Zhang 
40775024027SHong Zhang     ierr = KSPCreate(subcomm,&red->ksp);CHKERRQ(ierr);
40875024027SHong Zhang     ierr = KSPSetErrorIfNotConverged(red->ksp,pc->erroriffailure);CHKERRQ(ierr);
40975024027SHong Zhang     ierr = PetscObjectIncrementTabLevel((PetscObject)red->ksp,(PetscObject)pc,1);CHKERRQ(ierr);
41075024027SHong Zhang     ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)red->ksp);CHKERRQ(ierr);
41175024027SHong Zhang     ierr = KSPSetType(red->ksp,KSPPREONLY);CHKERRQ(ierr);
41275024027SHong Zhang     ierr = KSPGetPC(red->ksp,&red->pc);CHKERRQ(ierr);
41375024027SHong Zhang     ierr = PCSetType(red->pc,PCLU);CHKERRQ(ierr);
414753b7fb9SBarry Smith     if (red->shifttypeset) {
415753b7fb9SBarry Smith       ierr = PCFactorSetShiftType(red->pc,red->shifttype);CHKERRQ(ierr);
416753b7fb9SBarry Smith       red->shifttypeset = PETSC_FALSE;
417753b7fb9SBarry Smith     }
41875024027SHong Zhang     ierr = KSPSetOptionsPrefix(red->ksp,prefix);CHKERRQ(ierr);
41975024027SHong Zhang     ierr = KSPAppendOptionsPrefix(red->ksp,"redundant_");CHKERRQ(ierr);
42075024027SHong Zhang   }
42183ab6a24SBarry Smith   *innerksp = red->ksp;
4224b9ad928SBarry Smith   PetscFunctionReturn(0);
4234b9ad928SBarry Smith }
4244b9ad928SBarry Smith 
4254b9ad928SBarry Smith /*@
42683ab6a24SBarry Smith    PCRedundantGetKSP - Gets the less parallel KSP created by the redundant PC.
4274b9ad928SBarry Smith 
4284b9ad928SBarry Smith    Not Collective
4294b9ad928SBarry Smith 
4304b9ad928SBarry Smith    Input Parameter:
4314b9ad928SBarry Smith .  pc - the preconditioner context
4324b9ad928SBarry Smith 
4334b9ad928SBarry Smith    Output Parameter:
43483ab6a24SBarry Smith .  innerksp - the KSP on the smaller set of processes
4354b9ad928SBarry Smith 
4364b9ad928SBarry Smith    Level: advanced
4374b9ad928SBarry Smith 
4384b9ad928SBarry Smith .keywords: PC, redundant solve
4394b9ad928SBarry Smith @*/
44083ab6a24SBarry Smith PetscErrorCode PCRedundantGetKSP(PC pc,KSP *innerksp)
4414b9ad928SBarry Smith {
4424ac538c5SBarry Smith   PetscErrorCode ierr;
4434b9ad928SBarry Smith 
4444b9ad928SBarry Smith   PetscFunctionBegin;
4450700a824SBarry Smith   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
44683ab6a24SBarry Smith   PetscValidPointer(innerksp,2);
447753b7fb9SBarry Smith   ierr = PetscUseMethod(pc,"PCRedundantGetKSP_C",(PC,KSP*),(pc,innerksp));CHKERRQ(ierr);
4484b9ad928SBarry Smith   PetscFunctionReturn(0);
4494b9ad928SBarry Smith }
4504b9ad928SBarry Smith 
451f7a08781SBarry Smith static PetscErrorCode PCRedundantGetOperators_Redundant(PC pc,Mat *mat,Mat *pmat)
4524b9ad928SBarry Smith {
4534b9ad928SBarry Smith   PC_Redundant *red = (PC_Redundant*)pc->data;
4544b9ad928SBarry Smith 
4554b9ad928SBarry Smith   PetscFunctionBegin;
456b3804887SHong Zhang   if (mat)  *mat  = red->pmats;
457b3804887SHong Zhang   if (pmat) *pmat = red->pmats;
4584b9ad928SBarry Smith   PetscFunctionReturn(0);
4594b9ad928SBarry Smith }
4604b9ad928SBarry Smith 
4614b9ad928SBarry Smith /*@
4624b9ad928SBarry Smith    PCRedundantGetOperators - gets the sequential matrix and preconditioner matrix
4634b9ad928SBarry Smith 
4644b9ad928SBarry Smith    Not Collective
4654b9ad928SBarry Smith 
4664b9ad928SBarry Smith    Input Parameter:
4674b9ad928SBarry Smith .  pc - the preconditioner context
4684b9ad928SBarry Smith 
4694b9ad928SBarry Smith    Output Parameters:
4704b9ad928SBarry Smith +  mat - the matrix
4714b9ad928SBarry Smith -  pmat - the (possibly different) preconditioner matrix
4724b9ad928SBarry Smith 
4734b9ad928SBarry Smith    Level: advanced
4744b9ad928SBarry Smith 
4754b9ad928SBarry Smith .keywords: PC, redundant solve
4764b9ad928SBarry Smith @*/
4777087cfbeSBarry Smith PetscErrorCode PCRedundantGetOperators(PC pc,Mat *mat,Mat *pmat)
4784b9ad928SBarry Smith {
4794ac538c5SBarry Smith   PetscErrorCode ierr;
4804b9ad928SBarry Smith 
4814b9ad928SBarry Smith   PetscFunctionBegin;
4820700a824SBarry Smith   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
4834482741eSBarry Smith   if (mat)  PetscValidPointer(mat,2);
4844482741eSBarry Smith   if (pmat) PetscValidPointer(pmat,3);
485163d334eSBarry Smith   ierr = PetscUseMethod(pc,"PCRedundantGetOperators_C",(PC,Mat*,Mat*),(pc,mat,pmat));CHKERRQ(ierr);
4864b9ad928SBarry Smith   PetscFunctionReturn(0);
4874b9ad928SBarry Smith }
4884b9ad928SBarry Smith 
4894b9ad928SBarry Smith /* -------------------------------------------------------------------------------------*/
49037a17b4dSBarry Smith /*MC
49183ab6a24SBarry Smith      PCREDUNDANT - Runs a KSP solver with preconditioner for the entire problem on subgroups of processors
49237a17b4dSBarry Smith 
49383ab6a24SBarry Smith      Options for the redundant preconditioners can be set with -redundant_pc_xxx for the redundant KSP with -redundant_ksp_xxx
49437a17b4dSBarry Smith 
49509391456SBarry Smith   Options Database:
4969b21d695SBarry Smith .  -pc_redundant_number <n> - number of redundant solves, for example if you are using 64 MPI processes and
4979b21d695SBarry Smith                               use an n of 4 there will be 4 parallel solves each on 16 = 64/4 processes.
49809391456SBarry Smith 
49937a17b4dSBarry Smith    Level: intermediate
50037a17b4dSBarry Smith 
50195452b02SPatrick Sanan    Notes:
50295452b02SPatrick Sanan     The default KSP is preonly and the default PC is LU.
50383ab6a24SBarry Smith 
504753b7fb9SBarry Smith    PCFactorSetShiftType() applied to this PC will convey they shift type into the inner PC if it is factorization based.
505753b7fb9SBarry Smith 
50695452b02SPatrick Sanan    Developer Notes:
50795452b02SPatrick Sanan     Note that PCSetInitialGuessNonzero()  is not used by this class but likely should be.
5089cfaa89bSBarry Smith 
50937a17b4dSBarry Smith .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PCRedundantSetScatter(),
51083ab6a24SBarry Smith            PCRedundantGetKSP(), PCRedundantGetOperators(), PCRedundantSetNumber()
51137a17b4dSBarry Smith M*/
51237a17b4dSBarry Smith 
5138cc058d9SJed Brown PETSC_EXTERN PetscErrorCode PCCreate_Redundant(PC pc)
5144b9ad928SBarry Smith {
515dfbe8321SBarry Smith   PetscErrorCode ierr;
5164b9ad928SBarry Smith   PC_Redundant   *red;
51769db28dcSHong Zhang   PetscMPIInt    size;
5183f457be1SHong Zhang 
5194b9ad928SBarry Smith   PetscFunctionBegin;
520b00a9115SJed Brown   ierr = PetscNewLog(pc,&red);CHKERRQ(ierr);
521ce94432eSBarry Smith   ierr = MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);CHKERRQ(ierr);
5222fa5cd67SKarl Rupp 
52369db28dcSHong Zhang   red->nsubcomm       = size;
5244b9ad928SBarry Smith   red->useparallelmat = PETSC_TRUE;
5251fbd8f88SHong Zhang   pc->data            = (void*)red;
5264b9ad928SBarry Smith 
5274b9ad928SBarry Smith   pc->ops->apply          = PCApply_Redundant;
528d88bfacbSStefano Zampini   pc->ops->applytranspose = PCApplyTranspose_Redundant;
5294b9ad928SBarry Smith   pc->ops->setup          = PCSetUp_Redundant;
5304b9ad928SBarry Smith   pc->ops->destroy        = PCDestroy_Redundant;
5311ea5a559SBarry Smith   pc->ops->reset          = PCReset_Redundant;
5324b9ad928SBarry Smith   pc->ops->setfromoptions = PCSetFromOptions_Redundant;
5334b9ad928SBarry Smith   pc->ops->view           = PCView_Redundant;
5342fa5cd67SKarl Rupp 
535bdf89e91SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCRedundantSetScatter_C",PCRedundantSetScatter_Redundant);CHKERRQ(ierr);
536bdf89e91SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCRedundantSetNumber_C",PCRedundantSetNumber_Redundant);CHKERRQ(ierr);
537bdf89e91SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCRedundantGetKSP_C",PCRedundantGetKSP_Redundant);CHKERRQ(ierr);
538bdf89e91SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCRedundantGetOperators_C",PCRedundantGetOperators_Redundant);CHKERRQ(ierr);
539753b7fb9SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftType_C",PCFactorSetShiftType_Redundant);CHKERRQ(ierr);
5404b9ad928SBarry Smith   PetscFunctionReturn(0);
5414b9ad928SBarry Smith }
542b2573a8aSBarry Smith 
543