xref: /petsc/src/ksp/pc/impls/redundant/redundant.c (revision 2e956fe4fc852fabc23b437482e1fb7b77fddb0d)
1dba47a55SKris Buschelman 
24b9ad928SBarry Smith /*
33f457be1SHong Zhang   This file defines a "solve the problem redundantly on each subgroup of processor" preconditioner.
44b9ad928SBarry Smith */
5af0996ceSBarry Smith #include <petsc/private/pcimpl.h>
6c6db04a5SJed Brown #include <petscksp.h>           /*I "petscksp.h" I*/
74b9ad928SBarry Smith 
84b9ad928SBarry Smith typedef struct {
93e065800SHong Zhang   KSP                ksp;
104b9ad928SBarry Smith   PC                 pc;                   /* actual preconditioner used on each processor */
11ce94432eSBarry Smith   Vec                xsub,ysub;            /* vectors of a subcommunicator to hold parallel vectors of PetscObjectComm((PetscObject)pc) */
123f457be1SHong Zhang   Vec                xdup,ydup;            /* parallel vector that congregates xsub or ysub facilitating vector scattering */
13b3804887SHong Zhang   Mat                pmats;                /* matrix and optional preconditioner matrix belong to a subcommunicator */
143f457be1SHong Zhang   VecScatter         scatterin,scatterout; /* scatter used to move all values to each processor group (subcommunicator) */
15ace3abfcSBarry Smith   PetscBool          useparallelmat;
16c540e29cSHong Zhang   PetscSubcomm       psubcomm;
171fbd8f88SHong Zhang   PetscInt           nsubcomm;             /* num of data structure PetscSubcomm */
18753b7fb9SBarry Smith   PetscBool          shifttypeset;
19753b7fb9SBarry Smith   MatFactorShiftType shifttype;
204b9ad928SBarry Smith } PC_Redundant;
214b9ad928SBarry Smith 
22753b7fb9SBarry Smith PetscErrorCode  PCFactorSetShiftType_Redundant(PC pc,MatFactorShiftType shifttype)
23753b7fb9SBarry Smith {
24753b7fb9SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
25753b7fb9SBarry Smith 
26753b7fb9SBarry Smith   PetscFunctionBegin;
27753b7fb9SBarry Smith   if (red->ksp) {
28753b7fb9SBarry Smith     PC pc;
299566063dSJacob Faibussowitsch     PetscCall(KSPGetPC(red->ksp,&pc));
309566063dSJacob Faibussowitsch     PetscCall(PCFactorSetShiftType(pc,shifttype));
31753b7fb9SBarry Smith   } else {
32753b7fb9SBarry Smith     red->shifttypeset = PETSC_TRUE;
33753b7fb9SBarry Smith     red->shifttype    = shifttype;
34753b7fb9SBarry Smith   }
35753b7fb9SBarry Smith   PetscFunctionReturn(0);
36753b7fb9SBarry Smith }
37753b7fb9SBarry Smith 
386849ba73SBarry Smith static PetscErrorCode PCView_Redundant(PC pc,PetscViewer viewer)
394b9ad928SBarry Smith {
404b9ad928SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
41ace3abfcSBarry Smith   PetscBool      iascii,isstring;
4203ccd0b4SBarry Smith   PetscViewer    subviewer;
434b9ad928SBarry Smith 
444b9ad928SBarry Smith   PetscFunctionBegin;
459566063dSJacob Faibussowitsch   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii));
469566063dSJacob Faibussowitsch   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring));
4732077d6dSBarry Smith   if (iascii) {
4803ccd0b4SBarry Smith     if (!red->psubcomm) {
499566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(viewer,"  Not yet setup\n"));
5003ccd0b4SBarry Smith     } else {
5163a3b9bcSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(viewer,"  First (color=0) of %" PetscInt_FMT " PCs follows\n",red->nsubcomm));
529566063dSJacob Faibussowitsch       PetscCall(PetscViewerGetSubViewer(viewer,((PetscObject)red->pc)->comm,&subviewer));
53f5dd71faSBarry Smith       if (!red->psubcomm->color) { /* only view first redundant pc */
549566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPushTab(subviewer));
559566063dSJacob Faibussowitsch         PetscCall(KSPView(red->ksp,subviewer));
569566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPopTab(subviewer));
574b9ad928SBarry Smith       }
589566063dSJacob Faibussowitsch       PetscCall(PetscViewerRestoreSubViewer(viewer,((PetscObject)red->pc)->comm,&subviewer));
594b9ad928SBarry Smith     }
6003ccd0b4SBarry Smith   } else if (isstring) {
619566063dSJacob Faibussowitsch     PetscCall(PetscViewerStringSPrintf(viewer," Redundant solver preconditioner"));
624b9ad928SBarry Smith   }
634b9ad928SBarry Smith   PetscFunctionReturn(0);
644b9ad928SBarry Smith }
654b9ad928SBarry Smith 
6619b3b6edSHong Zhang #include <../src/mat/impls/aij/mpi/mpiaij.h>
676849ba73SBarry Smith static PetscErrorCode PCSetUp_Redundant(PC pc)
684b9ad928SBarry Smith {
694b9ad928SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
701b81debcSHong Zhang   PetscInt       mstart,mend,mlocal,M;
7113f74950SBarry Smith   PetscMPIInt    size;
72ce94432eSBarry Smith   MPI_Comm       comm,subcomm;
73ddc54837SHong Zhang   Vec            x;
743f457be1SHong Zhang 
754b9ad928SBarry Smith   PetscFunctionBegin;
769566063dSJacob Faibussowitsch   PetscCall(PetscObjectGetComm((PetscObject)pc,&comm));
77ddc54837SHong Zhang 
78ddc54837SHong Zhang   /* if pmatrix set by user is sequential then we do not need to gather the parallel matrix */
799566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_size(comm,&size));
80ddc54837SHong Zhang   if (size == 1) red->useparallelmat = PETSC_FALSE;
811fbd8f88SHong Zhang 
824b9ad928SBarry Smith   if (!pc->setupcalled) {
831b81debcSHong Zhang     PetscInt mloc_sub;
8475024027SHong Zhang     if (!red->psubcomm) { /* create red->psubcomm, new ksp and pc over subcomm */
8575024027SHong Zhang       KSP ksp;
869566063dSJacob Faibussowitsch       PetscCall(PCRedundantGetKSP(pc,&ksp));
871b81debcSHong Zhang     }
8875024027SHong Zhang     subcomm = PetscSubcommChild(red->psubcomm);
891fbd8f88SHong Zhang 
901b81debcSHong Zhang     if (red->useparallelmat) {
911b81debcSHong Zhang       /* grab the parallel matrix and put it into processors of a subcomminicator */
929566063dSJacob Faibussowitsch       PetscCall(MatCreateRedundantMatrix(pc->pmat,red->psubcomm->n,subcomm,MAT_INITIAL_MATRIX,&red->pmats));
93b85f2e9bSHong Zhang 
949566063dSJacob Faibussowitsch       PetscCallMPI(MPI_Comm_size(subcomm,&size));
95b85f2e9bSHong Zhang       if (size > 1) {
9608cecb0aSPierre Jolivet         PetscBool foundpack,issbaij;
979566063dSJacob Faibussowitsch         PetscCall(PetscObjectTypeCompare((PetscObject)red->pmats,MATMPISBAIJ,&issbaij));
9808cecb0aSPierre Jolivet         if (!issbaij) {
999566063dSJacob Faibussowitsch           PetscCall(MatGetFactorAvailable(red->pmats,NULL,MAT_FACTOR_LU,&foundpack));
10008cecb0aSPierre Jolivet         } else {
1019566063dSJacob Faibussowitsch           PetscCall(MatGetFactorAvailable(red->pmats,NULL,MAT_FACTOR_CHOLESKY,&foundpack));
10208cecb0aSPierre Jolivet         }
103b85f2e9bSHong Zhang         if (!foundpack) { /* reset default ksp and pc */
1049566063dSJacob Faibussowitsch           PetscCall(KSPSetType(red->ksp,KSPGMRES));
1059566063dSJacob Faibussowitsch           PetscCall(PCSetType(red->pc,PCBJACOBI));
106c1619fb6SBarry Smith         } else {
1079566063dSJacob Faibussowitsch           PetscCall(PCFactorSetMatSolverType(red->pc,NULL));
108b85f2e9bSHong Zhang         }
109b85f2e9bSHong Zhang       }
110b85f2e9bSHong Zhang 
1119566063dSJacob Faibussowitsch       PetscCall(KSPSetOperators(red->ksp,red->pmats,red->pmats));
1124b9ad928SBarry Smith 
1131b81debcSHong Zhang       /* get working vectors xsub and ysub */
1149566063dSJacob Faibussowitsch       PetscCall(MatCreateVecs(red->pmats,&red->xsub,&red->ysub));
1152fa5cd67SKarl Rupp 
1168b96b0d2SHong Zhang       /* create working vectors xdup and ydup.
1178b96b0d2SHong Zhang        xdup concatenates all xsub's contigously to form a mpi vector over dupcomm  (see PetscSubcommCreate_interlaced())
1188b96b0d2SHong Zhang        ydup concatenates all ysub and has empty local arrays because ysub's arrays will be place into it.
119ce94432eSBarry Smith        Note: we use communicator dupcomm, not PetscObjectComm((PetscObject)pc)! */
1209566063dSJacob Faibussowitsch       PetscCall(MatGetLocalSize(red->pmats,&mloc_sub,NULL));
1219566063dSJacob Faibussowitsch       PetscCall(VecCreateMPI(PetscSubcommContiguousParent(red->psubcomm),mloc_sub,PETSC_DECIDE,&red->xdup));
1229566063dSJacob Faibussowitsch       PetscCall(VecCreateMPIWithArray(PetscSubcommContiguousParent(red->psubcomm),1,mloc_sub,PETSC_DECIDE,NULL,&red->ydup));
1233f457be1SHong Zhang 
124f68be91cSHong Zhang       /* create vecscatters */
125f68be91cSHong Zhang       if (!red->scatterin) { /* efficiency of scatterin is independent from psubcomm_type! */
1263f457be1SHong Zhang         IS       is1,is2;
1273f457be1SHong Zhang         PetscInt *idx1,*idx2,i,j,k;
12845fc02eaSBarry Smith 
1299566063dSJacob Faibussowitsch         PetscCall(MatCreateVecs(pc->pmat,&x,NULL));
1309566063dSJacob Faibussowitsch         PetscCall(VecGetSize(x,&M));
1319566063dSJacob Faibussowitsch         PetscCall(VecGetOwnershipRange(x,&mstart,&mend));
1321b81debcSHong Zhang         mlocal = mend - mstart;
1339566063dSJacob Faibussowitsch         PetscCall(PetscMalloc2(red->psubcomm->n*mlocal,&idx1,red->psubcomm->n*mlocal,&idx2));
1343f457be1SHong Zhang         j    = 0;
1351fbd8f88SHong Zhang         for (k=0; k<red->psubcomm->n; k++) {
1363f457be1SHong Zhang           for (i=mstart; i<mend; i++) {
1373f457be1SHong Zhang             idx1[j]   = i;
138ddc54837SHong Zhang             idx2[j++] = i + M*k;
1393f457be1SHong Zhang           }
1403f457be1SHong Zhang         }
1419566063dSJacob Faibussowitsch         PetscCall(ISCreateGeneral(comm,red->psubcomm->n*mlocal,idx1,PETSC_COPY_VALUES,&is1));
1429566063dSJacob Faibussowitsch         PetscCall(ISCreateGeneral(comm,red->psubcomm->n*mlocal,idx2,PETSC_COPY_VALUES,&is2));
1439566063dSJacob Faibussowitsch         PetscCall(VecScatterCreate(x,is1,red->xdup,is2,&red->scatterin));
1449566063dSJacob Faibussowitsch         PetscCall(ISDestroy(&is1));
1459566063dSJacob Faibussowitsch         PetscCall(ISDestroy(&is2));
1463f457be1SHong Zhang 
1476909748bSHong Zhang         /* Impl below is good for PETSC_SUBCOMM_INTERLACED (no inter-process communication) and PETSC_SUBCOMM_CONTIGUOUS (communication within subcomm) */
1489566063dSJacob Faibussowitsch         PetscCall(ISCreateStride(comm,mlocal,mstart+ red->psubcomm->color*M,1,&is1));
1499566063dSJacob Faibussowitsch         PetscCall(ISCreateStride(comm,mlocal,mstart,1,&is2));
1509566063dSJacob Faibussowitsch         PetscCall(VecScatterCreate(red->xdup,is1,x,is2,&red->scatterout));
1519566063dSJacob Faibussowitsch         PetscCall(ISDestroy(&is1));
1529566063dSJacob Faibussowitsch         PetscCall(ISDestroy(&is2));
1539566063dSJacob Faibussowitsch         PetscCall(PetscFree2(idx1,idx2));
1549566063dSJacob Faibussowitsch         PetscCall(VecDestroy(&x));
1551b81debcSHong Zhang       }
156ab661555SHong Zhang     } else { /* !red->useparallelmat */
1579566063dSJacob Faibussowitsch       PetscCall(KSPSetOperators(red->ksp,pc->mat,pc->pmat));
1581b81debcSHong Zhang     }
159ab661555SHong Zhang   } else { /* pc->setupcalled */
1604b9ad928SBarry Smith     if (red->useparallelmat) {
161ab661555SHong Zhang       MatReuse       reuse;
1621b81debcSHong Zhang       /* grab the parallel matrix and put it into processors of a subcomminicator */
1631b81debcSHong Zhang       /*--------------------------------------------------------------------------*/
164ab661555SHong Zhang       if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1654b9ad928SBarry Smith         /* destroy old matrices */
1669566063dSJacob Faibussowitsch         PetscCall(MatDestroy(&red->pmats));
167ab661555SHong Zhang         reuse = MAT_INITIAL_MATRIX;
1684b9ad928SBarry Smith       } else {
169ab661555SHong Zhang         reuse = MAT_REUSE_MATRIX;
170ab661555SHong Zhang       }
1719566063dSJacob Faibussowitsch       PetscCall(MatCreateRedundantMatrix(pc->pmat,red->psubcomm->n,PetscSubcommChild(red->psubcomm),reuse,&red->pmats));
1729566063dSJacob Faibussowitsch       PetscCall(KSPSetOperators(red->ksp,red->pmats,red->pmats));
173ab661555SHong Zhang     } else { /* !red->useparallelmat */
1749566063dSJacob Faibussowitsch       PetscCall(KSPSetOperators(red->ksp,pc->mat,pc->pmat));
1754b9ad928SBarry Smith     }
176ab661555SHong Zhang   }
1771b81debcSHong Zhang 
1780c24e6a1SHong Zhang   if (pc->setfromoptionscalled) {
1799566063dSJacob Faibussowitsch     PetscCall(KSPSetFromOptions(red->ksp));
1800c24e6a1SHong Zhang   }
1819566063dSJacob Faibussowitsch   PetscCall(KSPSetUp(red->ksp));
1824b9ad928SBarry Smith   PetscFunctionReturn(0);
1834b9ad928SBarry Smith }
1844b9ad928SBarry Smith 
1856849ba73SBarry Smith static PetscErrorCode PCApply_Redundant(PC pc,Vec x,Vec y)
1864b9ad928SBarry Smith {
1874b9ad928SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
1883f457be1SHong Zhang   PetscScalar    *array;
1894b9ad928SBarry Smith 
1904b9ad928SBarry Smith   PetscFunctionBegin;
191ddc54837SHong Zhang   if (!red->useparallelmat) {
1929566063dSJacob Faibussowitsch     PetscCall(KSPSolve(red->ksp,x,y));
1939566063dSJacob Faibussowitsch     PetscCall(KSPCheckSolve(red->ksp,pc,y));
194ddc54837SHong Zhang     PetscFunctionReturn(0);
195ddc54837SHong Zhang   }
196ddc54837SHong Zhang 
1973f457be1SHong Zhang   /* scatter x to xdup */
1989566063dSJacob Faibussowitsch   PetscCall(VecScatterBegin(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD));
1999566063dSJacob Faibussowitsch   PetscCall(VecScatterEnd(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD));
2003f457be1SHong Zhang 
2013f457be1SHong Zhang   /* place xdup's local array into xsub */
2029566063dSJacob Faibussowitsch   PetscCall(VecGetArray(red->xdup,&array));
2039566063dSJacob Faibussowitsch   PetscCall(VecPlaceArray(red->xsub,(const PetscScalar*)array));
2044b9ad928SBarry Smith 
2054b9ad928SBarry Smith   /* apply preconditioner on each processor */
2069566063dSJacob Faibussowitsch   PetscCall(KSPSolve(red->ksp,red->xsub,red->ysub));
2079566063dSJacob Faibussowitsch   PetscCall(KSPCheckSolve(red->ksp,pc,red->ysub));
2089566063dSJacob Faibussowitsch   PetscCall(VecResetArray(red->xsub));
2099566063dSJacob Faibussowitsch   PetscCall(VecRestoreArray(red->xdup,&array));
2104b9ad928SBarry Smith 
2113f457be1SHong Zhang   /* place ysub's local array into ydup */
2129566063dSJacob Faibussowitsch   PetscCall(VecGetArray(red->ysub,&array));
2139566063dSJacob Faibussowitsch   PetscCall(VecPlaceArray(red->ydup,(const PetscScalar*)array));
2143f457be1SHong Zhang 
2153f457be1SHong Zhang   /* scatter ydup to y */
2169566063dSJacob Faibussowitsch   PetscCall(VecScatterBegin(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD));
2179566063dSJacob Faibussowitsch   PetscCall(VecScatterEnd(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD));
2189566063dSJacob Faibussowitsch   PetscCall(VecResetArray(red->ydup));
2199566063dSJacob Faibussowitsch   PetscCall(VecRestoreArray(red->ysub,&array));
2204b9ad928SBarry Smith   PetscFunctionReturn(0);
2214b9ad928SBarry Smith }
2224b9ad928SBarry Smith 
223d88bfacbSStefano Zampini static PetscErrorCode PCApplyTranspose_Redundant(PC pc,Vec x,Vec y)
224d88bfacbSStefano Zampini {
225d88bfacbSStefano Zampini   PC_Redundant   *red = (PC_Redundant*)pc->data;
226d88bfacbSStefano Zampini   PetscScalar    *array;
227d88bfacbSStefano Zampini 
228d88bfacbSStefano Zampini   PetscFunctionBegin;
229d88bfacbSStefano Zampini   if (!red->useparallelmat) {
2309566063dSJacob Faibussowitsch     PetscCall(KSPSolveTranspose(red->ksp,x,y));
2319566063dSJacob Faibussowitsch     PetscCall(KSPCheckSolve(red->ksp,pc,y));
232d88bfacbSStefano Zampini     PetscFunctionReturn(0);
233d88bfacbSStefano Zampini   }
234d88bfacbSStefano Zampini 
235d88bfacbSStefano Zampini   /* scatter x to xdup */
2369566063dSJacob Faibussowitsch   PetscCall(VecScatterBegin(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD));
2379566063dSJacob Faibussowitsch   PetscCall(VecScatterEnd(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD));
238d88bfacbSStefano Zampini 
239d88bfacbSStefano Zampini   /* place xdup's local array into xsub */
2409566063dSJacob Faibussowitsch   PetscCall(VecGetArray(red->xdup,&array));
2419566063dSJacob Faibussowitsch   PetscCall(VecPlaceArray(red->xsub,(const PetscScalar*)array));
242d88bfacbSStefano Zampini 
243d88bfacbSStefano Zampini   /* apply preconditioner on each processor */
2449566063dSJacob Faibussowitsch   PetscCall(KSPSolveTranspose(red->ksp,red->xsub,red->ysub));
2459566063dSJacob Faibussowitsch   PetscCall(KSPCheckSolve(red->ksp,pc,red->ysub));
2469566063dSJacob Faibussowitsch   PetscCall(VecResetArray(red->xsub));
2479566063dSJacob Faibussowitsch   PetscCall(VecRestoreArray(red->xdup,&array));
248d88bfacbSStefano Zampini 
249d88bfacbSStefano Zampini   /* place ysub's local array into ydup */
2509566063dSJacob Faibussowitsch   PetscCall(VecGetArray(red->ysub,&array));
2519566063dSJacob Faibussowitsch   PetscCall(VecPlaceArray(red->ydup,(const PetscScalar*)array));
252d88bfacbSStefano Zampini 
253d88bfacbSStefano Zampini   /* scatter ydup to y */
2549566063dSJacob Faibussowitsch   PetscCall(VecScatterBegin(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD));
2559566063dSJacob Faibussowitsch   PetscCall(VecScatterEnd(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD));
2569566063dSJacob Faibussowitsch   PetscCall(VecResetArray(red->ydup));
2579566063dSJacob Faibussowitsch   PetscCall(VecRestoreArray(red->ysub,&array));
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;
2644b9ad928SBarry Smith 
2654b9ad928SBarry Smith   PetscFunctionBegin;
2661b81debcSHong Zhang   if (red->useparallelmat) {
2679566063dSJacob Faibussowitsch     PetscCall(VecScatterDestroy(&red->scatterin));
2689566063dSJacob Faibussowitsch     PetscCall(VecScatterDestroy(&red->scatterout));
2699566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&red->ysub));
2709566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&red->xsub));
2719566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&red->xdup));
2729566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&red->ydup));
2731b81debcSHong Zhang   }
2749566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&red->pmats));
2759566063dSJacob Faibussowitsch   PetscCall(KSPReset(red->ksp));
2761ea5a559SBarry Smith   PetscFunctionReturn(0);
2771ea5a559SBarry Smith }
2781ea5a559SBarry Smith 
2791ea5a559SBarry Smith static PetscErrorCode PCDestroy_Redundant(PC pc)
2801ea5a559SBarry Smith {
2811ea5a559SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
2821ea5a559SBarry Smith 
2831ea5a559SBarry Smith   PetscFunctionBegin;
2849566063dSJacob Faibussowitsch   PetscCall(PCReset_Redundant(pc));
2859566063dSJacob Faibussowitsch   PetscCall(KSPDestroy(&red->ksp));
2869566063dSJacob Faibussowitsch   PetscCall(PetscSubcommDestroy(&red->psubcomm));
287*2e956fe4SStefano Zampini   PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCRedundantSetScatter_C",NULL));
288*2e956fe4SStefano Zampini   PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCRedundantSetNumber_C",NULL));
289*2e956fe4SStefano Zampini   PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCRedundantGetKSP_C",NULL));
290*2e956fe4SStefano Zampini   PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCRedundantGetOperators_C",NULL));
291*2e956fe4SStefano Zampini   PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftType_C",NULL));
2929566063dSJacob Faibussowitsch   PetscCall(PetscFree(pc->data));
2934b9ad928SBarry Smith   PetscFunctionReturn(0);
2944b9ad928SBarry Smith }
2954b9ad928SBarry Smith 
2964416b707SBarry Smith static PetscErrorCode PCSetFromOptions_Redundant(PetscOptionItems *PetscOptionsObject,PC pc)
2974b9ad928SBarry Smith {
298a98ce0f4SHong Zhang   PC_Redundant   *red = (PC_Redundant*)pc->data;
299a98ce0f4SHong Zhang 
3004b9ad928SBarry Smith   PetscFunctionBegin;
301d0609cedSBarry Smith   PetscOptionsHeadBegin(PetscOptionsObject,"Redundant options");
3029566063dSJacob Faibussowitsch   PetscCall(PetscOptionsInt("-pc_redundant_number","Number of redundant pc","PCRedundantSetNumber",red->nsubcomm,&red->nsubcomm,NULL));
303d0609cedSBarry Smith   PetscOptionsHeadEnd();
3044b9ad928SBarry Smith   PetscFunctionReturn(0);
3054b9ad928SBarry Smith }
3064b9ad928SBarry Smith 
307f7a08781SBarry Smith static PetscErrorCode PCRedundantSetNumber_Redundant(PC pc,PetscInt nreds)
30809a6bc64SHong Zhang {
30909a6bc64SHong Zhang   PC_Redundant *red = (PC_Redundant*)pc->data;
31009a6bc64SHong Zhang 
31109a6bc64SHong Zhang   PetscFunctionBegin;
31209a6bc64SHong Zhang   red->nsubcomm = nreds;
31309a6bc64SHong Zhang   PetscFunctionReturn(0);
31409a6bc64SHong Zhang }
31509a6bc64SHong Zhang 
31609a6bc64SHong Zhang /*@
31709a6bc64SHong Zhang    PCRedundantSetNumber - Sets the number of redundant preconditioner contexts.
31809a6bc64SHong Zhang 
3193f9fe445SBarry Smith    Logically Collective on PC
32009a6bc64SHong Zhang 
32109a6bc64SHong Zhang    Input Parameters:
32209a6bc64SHong Zhang +  pc - the preconditioner context
3239b21d695SBarry Smith -  nredundant - number of redundant preconditioner contexts; for example if you are using 64 MPI processes and
3249b21d695SBarry Smith                               use an nredundant of 4 there will be 4 parallel solves each on 16 = 64/4 processes.
32509a6bc64SHong Zhang 
32609a6bc64SHong Zhang    Level: advanced
32709a6bc64SHong Zhang 
32809a6bc64SHong Zhang @*/
3297087cfbeSBarry Smith PetscErrorCode PCRedundantSetNumber(PC pc,PetscInt nredundant)
33009a6bc64SHong Zhang {
33109a6bc64SHong Zhang   PetscFunctionBegin;
3320700a824SBarry Smith   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
33363a3b9bcSJacob Faibussowitsch   PetscCheck(nredundant > 0,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "num of redundant pc %" PetscInt_FMT " must be positive",nredundant);
334cac4c232SBarry Smith   PetscTryMethod(pc,"PCRedundantSetNumber_C",(PC,PetscInt),(pc,nredundant));
33509a6bc64SHong Zhang   PetscFunctionReturn(0);
33609a6bc64SHong Zhang }
33709a6bc64SHong Zhang 
338f7a08781SBarry Smith static PetscErrorCode PCRedundantSetScatter_Redundant(PC pc,VecScatter in,VecScatter out)
3394b9ad928SBarry Smith {
3404b9ad928SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
3414b9ad928SBarry Smith 
3424b9ad928SBarry Smith   PetscFunctionBegin;
3439566063dSJacob Faibussowitsch   PetscCall(PetscObjectReference((PetscObject)in));
3449566063dSJacob Faibussowitsch   PetscCall(VecScatterDestroy(&red->scatterin));
3452fa5cd67SKarl Rupp 
346c3122656SLisandro Dalcin   red->scatterin  = in;
3472fa5cd67SKarl Rupp 
3489566063dSJacob Faibussowitsch   PetscCall(PetscObjectReference((PetscObject)out));
3499566063dSJacob Faibussowitsch   PetscCall(VecScatterDestroy(&red->scatterout));
350c3122656SLisandro Dalcin   red->scatterout = out;
3514b9ad928SBarry Smith   PetscFunctionReturn(0);
3524b9ad928SBarry Smith }
3534b9ad928SBarry Smith 
3544b9ad928SBarry Smith /*@
3554b9ad928SBarry Smith    PCRedundantSetScatter - Sets the scatter used to copy values into the
3564b9ad928SBarry Smith      redundant local solve and the scatter to move them back into the global
3574b9ad928SBarry Smith      vector.
3584b9ad928SBarry Smith 
3593f9fe445SBarry Smith    Logically Collective on PC
3604b9ad928SBarry Smith 
3614b9ad928SBarry Smith    Input Parameters:
3624b9ad928SBarry Smith +  pc - the preconditioner context
3634b9ad928SBarry Smith .  in - the scatter to move the values in
3644b9ad928SBarry Smith -  out - the scatter to move them out
3654b9ad928SBarry Smith 
3664b9ad928SBarry Smith    Level: advanced
3674b9ad928SBarry Smith 
3684b9ad928SBarry Smith @*/
3697087cfbeSBarry Smith PetscErrorCode PCRedundantSetScatter(PC pc,VecScatter in,VecScatter out)
3704b9ad928SBarry Smith {
3714b9ad928SBarry Smith   PetscFunctionBegin;
3720700a824SBarry Smith   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
37397929ea7SJunchao Zhang   PetscValidHeaderSpecific(in,PETSCSF_CLASSID,2);
37497929ea7SJunchao Zhang   PetscValidHeaderSpecific(out,PETSCSF_CLASSID,3);
375cac4c232SBarry Smith   PetscTryMethod(pc,"PCRedundantSetScatter_C",(PC,VecScatter,VecScatter),(pc,in,out));
3764b9ad928SBarry Smith   PetscFunctionReturn(0);
3774b9ad928SBarry Smith }
3784b9ad928SBarry Smith 
379f7a08781SBarry Smith static PetscErrorCode PCRedundantGetKSP_Redundant(PC pc,KSP *innerksp)
3804b9ad928SBarry Smith {
3814b9ad928SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
38275024027SHong Zhang   MPI_Comm       comm,subcomm;
38375024027SHong Zhang   const char     *prefix;
38408cecb0aSPierre Jolivet   PetscBool      issbaij;
3854b9ad928SBarry Smith 
3864b9ad928SBarry Smith   PetscFunctionBegin;
38775024027SHong Zhang   if (!red->psubcomm) {
3889566063dSJacob Faibussowitsch     PetscCall(PCGetOptionsPrefix(pc,&prefix));
389e5acf8a4SHong Zhang 
3909566063dSJacob Faibussowitsch     PetscCall(PetscObjectGetComm((PetscObject)pc,&comm));
3919566063dSJacob Faibussowitsch     PetscCall(PetscSubcommCreate(comm,&red->psubcomm));
3929566063dSJacob Faibussowitsch     PetscCall(PetscSubcommSetNumber(red->psubcomm,red->nsubcomm));
3939566063dSJacob Faibussowitsch     PetscCall(PetscSubcommSetType(red->psubcomm,PETSC_SUBCOMM_CONTIGUOUS));
394e5acf8a4SHong Zhang 
3959566063dSJacob Faibussowitsch     PetscCall(PetscSubcommSetOptionsPrefix(red->psubcomm,prefix));
3969566063dSJacob Faibussowitsch     PetscCall(PetscSubcommSetFromOptions(red->psubcomm));
3979566063dSJacob Faibussowitsch     PetscCall(PetscLogObjectMemory((PetscObject)pc,sizeof(PetscSubcomm)));
39875024027SHong Zhang 
39975024027SHong Zhang     /* create a new PC that processors in each subcomm have copy of */
40075024027SHong Zhang     subcomm = PetscSubcommChild(red->psubcomm);
40175024027SHong Zhang 
4029566063dSJacob Faibussowitsch     PetscCall(KSPCreate(subcomm,&red->ksp));
4039566063dSJacob Faibussowitsch     PetscCall(KSPSetErrorIfNotConverged(red->ksp,pc->erroriffailure));
4049566063dSJacob Faibussowitsch     PetscCall(PetscObjectIncrementTabLevel((PetscObject)red->ksp,(PetscObject)pc,1));
4059566063dSJacob Faibussowitsch     PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)red->ksp));
4069566063dSJacob Faibussowitsch     PetscCall(KSPSetType(red->ksp,KSPPREONLY));
4079566063dSJacob Faibussowitsch     PetscCall(KSPGetPC(red->ksp,&red->pc));
4089566063dSJacob Faibussowitsch     PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat,MATSEQSBAIJ,&issbaij));
40908cecb0aSPierre Jolivet     if (!issbaij) {
4109566063dSJacob Faibussowitsch       PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat,MATMPISBAIJ,&issbaij));
41108cecb0aSPierre Jolivet     }
41208cecb0aSPierre Jolivet     if (!issbaij) {
4139566063dSJacob Faibussowitsch       PetscCall(PCSetType(red->pc,PCLU));
41408cecb0aSPierre Jolivet     } else {
4159566063dSJacob Faibussowitsch       PetscCall(PCSetType(red->pc,PCCHOLESKY));
41608cecb0aSPierre Jolivet     }
417753b7fb9SBarry Smith     if (red->shifttypeset) {
4189566063dSJacob Faibussowitsch       PetscCall(PCFactorSetShiftType(red->pc,red->shifttype));
419753b7fb9SBarry Smith       red->shifttypeset = PETSC_FALSE;
420753b7fb9SBarry Smith     }
4219566063dSJacob Faibussowitsch     PetscCall(KSPSetOptionsPrefix(red->ksp,prefix));
4229566063dSJacob Faibussowitsch     PetscCall(KSPAppendOptionsPrefix(red->ksp,"redundant_"));
42375024027SHong Zhang   }
42483ab6a24SBarry Smith   *innerksp = red->ksp;
4254b9ad928SBarry Smith   PetscFunctionReturn(0);
4264b9ad928SBarry Smith }
4274b9ad928SBarry Smith 
4284b9ad928SBarry Smith /*@
42983ab6a24SBarry Smith    PCRedundantGetKSP - Gets the less parallel KSP created by the redundant PC.
4304b9ad928SBarry Smith 
4314b9ad928SBarry Smith    Not Collective
4324b9ad928SBarry Smith 
4334b9ad928SBarry Smith    Input Parameter:
4344b9ad928SBarry Smith .  pc - the preconditioner context
4354b9ad928SBarry Smith 
4364b9ad928SBarry Smith    Output Parameter:
43783ab6a24SBarry Smith .  innerksp - the KSP on the smaller set of processes
4384b9ad928SBarry Smith 
4394b9ad928SBarry Smith    Level: advanced
4404b9ad928SBarry Smith 
4414b9ad928SBarry Smith @*/
44283ab6a24SBarry Smith PetscErrorCode PCRedundantGetKSP(PC pc,KSP *innerksp)
4434b9ad928SBarry Smith {
4444b9ad928SBarry Smith   PetscFunctionBegin;
4450700a824SBarry Smith   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
44683ab6a24SBarry Smith   PetscValidPointer(innerksp,2);
447cac4c232SBarry Smith   PetscUseMethod(pc,"PCRedundantGetKSP_C",(PC,KSP*),(pc,innerksp));
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 @*/
4767087cfbeSBarry Smith PetscErrorCode PCRedundantGetOperators(PC pc,Mat *mat,Mat *pmat)
4774b9ad928SBarry Smith {
4784b9ad928SBarry Smith   PetscFunctionBegin;
4790700a824SBarry Smith   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
4804482741eSBarry Smith   if (mat)  PetscValidPointer(mat,2);
4814482741eSBarry Smith   if (pmat) PetscValidPointer(pmat,3);
482cac4c232SBarry Smith   PetscUseMethod(pc,"PCRedundantGetOperators_C",(PC,Mat*,Mat*),(pc,mat,pmat));
4834b9ad928SBarry Smith   PetscFunctionReturn(0);
4844b9ad928SBarry Smith }
4854b9ad928SBarry Smith 
4864b9ad928SBarry Smith /* -------------------------------------------------------------------------------------*/
48737a17b4dSBarry Smith /*MC
48883ab6a24SBarry Smith      PCREDUNDANT - Runs a KSP solver with preconditioner for the entire problem on subgroups of processors
48937a17b4dSBarry Smith 
49083ab6a24SBarry Smith      Options for the redundant preconditioners can be set with -redundant_pc_xxx for the redundant KSP with -redundant_ksp_xxx
49137a17b4dSBarry Smith 
49209391456SBarry Smith   Options Database:
4939b21d695SBarry Smith .  -pc_redundant_number <n> - number of redundant solves, for example if you are using 64 MPI processes and
4949b21d695SBarry Smith                               use an n of 4 there will be 4 parallel solves each on 16 = 64/4 processes.
49509391456SBarry Smith 
49637a17b4dSBarry Smith    Level: intermediate
49737a17b4dSBarry Smith 
49895452b02SPatrick Sanan    Notes:
49908cecb0aSPierre Jolivet     The default KSP is preonly and the default PC is LU or CHOLESKY if Pmat is of type MATSBAIJ.
50083ab6a24SBarry Smith 
501753b7fb9SBarry Smith    PCFactorSetShiftType() applied to this PC will convey they shift type into the inner PC if it is factorization based.
502753b7fb9SBarry Smith 
50395452b02SPatrick Sanan    Developer Notes:
50495452b02SPatrick Sanan     Note that PCSetInitialGuessNonzero()  is not used by this class but likely should be.
5059cfaa89bSBarry Smith 
506db781477SPatrick Sanan .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PCRedundantSetScatter()`,
507db781477SPatrick Sanan           `PCRedundantGetKSP()`, `PCRedundantGetOperators()`, `PCRedundantSetNumber()`
50837a17b4dSBarry Smith M*/
50937a17b4dSBarry Smith 
5108cc058d9SJed Brown PETSC_EXTERN PetscErrorCode PCCreate_Redundant(PC pc)
5114b9ad928SBarry Smith {
5124b9ad928SBarry Smith   PC_Redundant   *red;
51369db28dcSHong Zhang   PetscMPIInt    size;
5143f457be1SHong Zhang 
5154b9ad928SBarry Smith   PetscFunctionBegin;
5169566063dSJacob Faibussowitsch   PetscCall(PetscNewLog(pc,&red));
5179566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size));
5182fa5cd67SKarl Rupp 
51969db28dcSHong Zhang   red->nsubcomm       = size;
5204b9ad928SBarry Smith   red->useparallelmat = PETSC_TRUE;
5211fbd8f88SHong Zhang   pc->data            = (void*)red;
5224b9ad928SBarry Smith 
5234b9ad928SBarry Smith   pc->ops->apply          = PCApply_Redundant;
524d88bfacbSStefano Zampini   pc->ops->applytranspose = PCApplyTranspose_Redundant;
5254b9ad928SBarry Smith   pc->ops->setup          = PCSetUp_Redundant;
5264b9ad928SBarry Smith   pc->ops->destroy        = PCDestroy_Redundant;
5271ea5a559SBarry Smith   pc->ops->reset          = PCReset_Redundant;
5284b9ad928SBarry Smith   pc->ops->setfromoptions = PCSetFromOptions_Redundant;
5294b9ad928SBarry Smith   pc->ops->view           = PCView_Redundant;
5302fa5cd67SKarl Rupp 
5319566063dSJacob Faibussowitsch   PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCRedundantSetScatter_C",PCRedundantSetScatter_Redundant));
5329566063dSJacob Faibussowitsch   PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCRedundantSetNumber_C",PCRedundantSetNumber_Redundant));
5339566063dSJacob Faibussowitsch   PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCRedundantGetKSP_C",PCRedundantGetKSP_Redundant));
5349566063dSJacob Faibussowitsch   PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCRedundantGetOperators_C",PCRedundantGetOperators_Redundant));
5359566063dSJacob Faibussowitsch   PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftType_C",PCFactorSetShiftType_Redundant));
5364b9ad928SBarry Smith   PetscFunctionReturn(0);
5374b9ad928SBarry Smith }
538