xref: /petsc/src/ksp/pc/impls/redundant/redundant.c (revision efd4aadf157bf1ba2d80c2be092fcf4247860003)
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