xref: /petsc/src/ksp/pc/impls/redundant/redundant.c (revision ddc54837807eb96c5bc9ef34be37eed03061befe)
1dba47a55SKris Buschelman 
24b9ad928SBarry Smith /*
33f457be1SHong Zhang   This file defines a "solve the problem redundantly on each subgroup of processor" preconditioner.
44b9ad928SBarry Smith */
507475bc1SBarry 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 */
184b9ad928SBarry Smith } PC_Redundant;
194b9ad928SBarry Smith 
204b9ad928SBarry Smith #undef __FUNCT__
214b9ad928SBarry Smith #define __FUNCT__ "PCView_Redundant"
226849ba73SBarry Smith static PetscErrorCode PCView_Redundant(PC pc,PetscViewer viewer)
234b9ad928SBarry Smith {
244b9ad928SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
25dfbe8321SBarry Smith   PetscErrorCode ierr;
26ace3abfcSBarry Smith   PetscBool      iascii,isstring;
2703ccd0b4SBarry Smith   PetscViewer    subviewer;
284b9ad928SBarry Smith 
294b9ad928SBarry Smith   PetscFunctionBegin;
30251f4c67SDmitry Karpeev   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
31251f4c67SDmitry Karpeev   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);CHKERRQ(ierr);
3232077d6dSBarry Smith   if (iascii) {
3303ccd0b4SBarry Smith     if (!red->psubcomm) {
3403ccd0b4SBarry Smith       ierr = PetscViewerASCIIPrintf(viewer,"  Redundant preconditioner: Not yet setup\n");CHKERRQ(ierr);
3503ccd0b4SBarry Smith     } else {
363e065800SHong Zhang       ierr = PetscViewerASCIIPrintf(viewer,"  Redundant preconditioner: First (color=0) of %D PCs follows\n",red->nsubcomm);CHKERRQ(ierr);
377adad957SLisandro Dalcin       ierr = PetscViewerGetSubcomm(viewer,((PetscObject)red->pc)->comm,&subviewer);CHKERRQ(ierr);
38f5dd71faSBarry Smith       if (!red->psubcomm->color) { /* only view first redundant pc */
394b9ad928SBarry Smith         ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
403e065800SHong Zhang         ierr = KSPView(red->ksp,subviewer);CHKERRQ(ierr);
414b9ad928SBarry Smith         ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
424b9ad928SBarry Smith       }
437adad957SLisandro Dalcin       ierr = PetscViewerRestoreSubcomm(viewer,((PetscObject)red->pc)->comm,&subviewer);CHKERRQ(ierr);
444b9ad928SBarry Smith     }
4503ccd0b4SBarry Smith   } else if (isstring) {
4603ccd0b4SBarry Smith     ierr = PetscViewerStringSPrintf(viewer," Redundant solver preconditioner");CHKERRQ(ierr);
474b9ad928SBarry Smith   }
484b9ad928SBarry Smith   PetscFunctionReturn(0);
494b9ad928SBarry Smith }
504b9ad928SBarry Smith 
51d3b23db5SHong Zhang extern PetscErrorCode MatGetRedundantMatrix_MPIAIJ_interlaced(Mat,PetscInt,PetscSubcomm,MatReuse,Mat*); /* rm later! */
52d3b23db5SHong Zhang 
534b9ad928SBarry Smith #undef __FUNCT__
544b9ad928SBarry Smith #define __FUNCT__ "PCSetUp_Redundant"
556849ba73SBarry Smith static PetscErrorCode PCSetUp_Redundant(PC pc)
564b9ad928SBarry Smith {
574b9ad928SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
58dfbe8321SBarry Smith   PetscErrorCode ierr;
59*ddc54837SHong Zhang   PetscInt       mstart,mend,mlocal,M,mlocal_sub,rstart_sub,rend_sub,mloc_sub;
6013f74950SBarry Smith   PetscMPIInt    size;
614b9ad928SBarry Smith   MatReuse       reuse = MAT_INITIAL_MATRIX;
624b9ad928SBarry Smith   MatStructure   str   = DIFFERENT_NONZERO_PATTERN;
63ce94432eSBarry Smith   MPI_Comm       comm,subcomm;
64*ddc54837SHong Zhang   Vec            x;
653f457be1SHong Zhang   PetscMPIInt    subsize,subrank;
661fbd8f88SHong Zhang   const char     *prefix;
67b862ddfaSBarry Smith   const PetscInt *range;
683f457be1SHong Zhang 
694b9ad928SBarry Smith   PetscFunctionBegin;
70ce94432eSBarry Smith   ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
71*ddc54837SHong Zhang   ierr = MatGetVecs(pc->pmat,&x,0);CHKERRQ(ierr);
72*ddc54837SHong Zhang   ierr = VecGetSize(x,&M);CHKERRQ(ierr);
73*ddc54837SHong Zhang 
74*ddc54837SHong Zhang   /* if pmatrix set by user is sequential then we do not need to gather the parallel matrix */
75*ddc54837SHong Zhang   ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
76*ddc54837SHong Zhang   if (size == 1) red->useparallelmat = PETSC_FALSE;
771fbd8f88SHong Zhang 
784b9ad928SBarry Smith   if (!pc->setupcalled) {
795f06b7aaSBarry Smith     if (!red->psubcomm) {
80d8a68f86SHong Zhang       ierr = PetscSubcommCreate(comm,&red->psubcomm);CHKERRQ(ierr);
81d8a68f86SHong Zhang       ierr = PetscSubcommSetNumber(red->psubcomm,red->nsubcomm);CHKERRQ(ierr);
82d8a68f86SHong Zhang       ierr = PetscSubcommSetType(red->psubcomm,PETSC_SUBCOMM_INTERLACED);CHKERRQ(ierr);
831fbd8f88SHong Zhang       ierr = PetscLogObjectMemory(pc,sizeof(PetscSubcomm));CHKERRQ(ierr);
841fbd8f88SHong Zhang 
851fbd8f88SHong Zhang       /* create a new PC that processors in each subcomm have copy of */
860d7810c8SBarry Smith       subcomm = red->psubcomm->comm;
872fa5cd67SKarl Rupp 
885f06b7aaSBarry Smith       ierr = KSPCreate(subcomm,&red->ksp);CHKERRQ(ierr);
895f06b7aaSBarry Smith       ierr = PetscObjectIncrementTabLevel((PetscObject)red->ksp,(PetscObject)pc,1);CHKERRQ(ierr);
905f06b7aaSBarry Smith       ierr = PetscLogObjectParent(pc,red->ksp);CHKERRQ(ierr);
915f06b7aaSBarry Smith       ierr = KSPSetType(red->ksp,KSPPREONLY);CHKERRQ(ierr);
925f06b7aaSBarry Smith       ierr = KSPGetPC(red->ksp,&red->pc);CHKERRQ(ierr);
93cf52b8b1SHong Zhang       ierr = PCSetType(red->pc,PCLU);CHKERRQ(ierr);
94cf52b8b1SHong Zhang 
951fbd8f88SHong Zhang       ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr);
965f06b7aaSBarry Smith       ierr = KSPSetOptionsPrefix(red->ksp,prefix);CHKERRQ(ierr);
975f06b7aaSBarry Smith       ierr = KSPAppendOptionsPrefix(red->ksp,"redundant_");CHKERRQ(ierr);
982fa5cd67SKarl Rupp     } else subcomm = red->psubcomm->comm;
991fbd8f88SHong Zhang 
1003f457be1SHong Zhang     /* create working vectors xsub/ysub and xdup/ydup */
101*ddc54837SHong Zhang     ierr = VecGetLocalSize(x,&mlocal);CHKERRQ(ierr);
102*ddc54837SHong Zhang     ierr = VecGetOwnershipRange(x,&mstart,&mend);CHKERRQ(ierr);
1034b9ad928SBarry Smith 
1043f457be1SHong Zhang     /* get local size of xsub/ysub */
1051fbd8f88SHong Zhang     ierr = MPI_Comm_size(subcomm,&subsize);CHKERRQ(ierr);
1061fbd8f88SHong Zhang     ierr = MPI_Comm_rank(subcomm,&subrank);CHKERRQ(ierr);
107b862ddfaSBarry Smith     ierr = MatGetOwnershipRanges(pc->pmat,&range);CHKERRQ(ierr);
108b862ddfaSBarry Smith     rstart_sub = range[red->psubcomm->n*subrank]; /* rstart in xsub/ysub */
1092fa5cd67SKarl Rupp     if (subrank+1 < subsize) rend_sub = range[red->psubcomm->n*(subrank+1)];
110*ddc54837SHong Zhang     else rend_sub = M;
1112fa5cd67SKarl Rupp 
1123f457be1SHong Zhang     mloc_sub = rend_sub - rstart_sub;
1131fbd8f88SHong Zhang     ierr     = VecCreateMPI(subcomm,mloc_sub,PETSC_DECIDE,&red->ysub);CHKERRQ(ierr);
1143f457be1SHong Zhang     /* create xsub with empty local arrays, because xdup's arrays will be placed into it */
1150298fd71SBarry Smith     ierr = VecCreateMPIWithArray(subcomm,1,mloc_sub,PETSC_DECIDE,NULL,&red->xsub);CHKERRQ(ierr);
1163f457be1SHong Zhang 
1173f457be1SHong Zhang     /* create xdup and ydup. ydup has empty local arrays because ysub's arrays will be place into it.
118ce94432eSBarry Smith        Note: we use communicator dupcomm, not PetscObjectComm((PetscObject)pc)! */
1191fbd8f88SHong Zhang     ierr = VecCreateMPI(red->psubcomm->dupparent,mloc_sub,PETSC_DECIDE,&red->xdup);CHKERRQ(ierr);
1200298fd71SBarry Smith     ierr = VecCreateMPIWithArray(red->psubcomm->dupparent,1,mloc_sub,PETSC_DECIDE,NULL,&red->ydup);CHKERRQ(ierr);
1213f457be1SHong Zhang 
1223f457be1SHong Zhang     /* create vec scatters */
1233f457be1SHong Zhang     if (!red->scatterin) {
1243f457be1SHong Zhang       IS       is1,is2;
1253f457be1SHong Zhang       PetscInt *idx1,*idx2,i,j,k;
12645fc02eaSBarry Smith 
1271d79065fSBarry Smith       ierr = PetscMalloc2(red->psubcomm->n*mlocal,PetscInt,&idx1,red->psubcomm->n*mlocal,PetscInt,&idx2);CHKERRQ(ierr);
1283f457be1SHong Zhang       j    = 0;
1291fbd8f88SHong Zhang       for (k=0; k<red->psubcomm->n; k++) {
1303f457be1SHong Zhang         for (i=mstart; i<mend; i++) {
1313f457be1SHong Zhang           idx1[j]   = i;
132*ddc54837SHong Zhang           idx2[j++] = i + M*k;
1333f457be1SHong Zhang         }
1343f457be1SHong Zhang       }
13570b3c8c7SBarry Smith       ierr = ISCreateGeneral(comm,red->psubcomm->n*mlocal,idx1,PETSC_COPY_VALUES,&is1);CHKERRQ(ierr);
13670b3c8c7SBarry Smith       ierr = ISCreateGeneral(comm,red->psubcomm->n*mlocal,idx2,PETSC_COPY_VALUES,&is2);CHKERRQ(ierr);
137*ddc54837SHong Zhang       ierr = VecScatterCreate(x,is1,red->xdup,is2,&red->scatterin);CHKERRQ(ierr);
138fcfd50ebSBarry Smith       ierr = ISDestroy(&is1);CHKERRQ(ierr);
139fcfd50ebSBarry Smith       ierr = ISDestroy(&is2);CHKERRQ(ierr);
1403f457be1SHong Zhang 
141*ddc54837SHong Zhang       ierr = ISCreateStride(comm,mlocal,mstart+ red->psubcomm->color*M,1,&is1);CHKERRQ(ierr);
1423f457be1SHong Zhang       ierr = ISCreateStride(comm,mlocal,mstart,1,&is2);CHKERRQ(ierr);
143*ddc54837SHong Zhang       ierr = VecScatterCreate(red->xdup,is1,x,is2,&red->scatterout);CHKERRQ(ierr);
144fcfd50ebSBarry Smith       ierr = ISDestroy(&is1);CHKERRQ(ierr);
145fcfd50ebSBarry Smith       ierr = ISDestroy(&is2);CHKERRQ(ierr);
1461d79065fSBarry Smith       ierr = PetscFree2(idx1,idx2);CHKERRQ(ierr);
1474b9ad928SBarry Smith     }
1484b9ad928SBarry Smith   }
149*ddc54837SHong Zhang   ierr = VecDestroy(&x);CHKERRQ(ierr);
1504b9ad928SBarry Smith 
1514b9ad928SBarry Smith   if (red->useparallelmat) {
1524b9ad928SBarry Smith     if (pc->setupcalled == 1 && pc->flag == DIFFERENT_NONZERO_PATTERN) {
1534b9ad928SBarry Smith       /* destroy old matrices */
1546bf464f9SBarry Smith       ierr = MatDestroy(&red->pmats);CHKERRQ(ierr);
1554b9ad928SBarry Smith     } else if (pc->setupcalled == 1) {
1564b9ad928SBarry Smith       reuse = MAT_REUSE_MATRIX;
1574b9ad928SBarry Smith       str   = SAME_NONZERO_PATTERN;
1584b9ad928SBarry Smith     }
1594b9ad928SBarry Smith 
1603f457be1SHong Zhang     /* grab the parallel matrix and put it into processors of a subcomminicator */
161f664ae05SHong Zhang     /*--------------------------------------------------------------------------*/
162f664ae05SHong Zhang     ierr = VecGetLocalSize(red->ysub,&mlocal_sub);CHKERRQ(ierr);
163d3b23db5SHong Zhang     ierr = MatGetRedundantMatrix_MPIAIJ_interlaced(pc->pmat,red->psubcomm->n,red->psubcomm,reuse,&red->pmats);CHKERRQ(ierr);
1643f457be1SHong Zhang     /* tell PC of the subcommunicator its operators */
16590f1c854SHong Zhang     ierr = KSPSetOperators(red->ksp,red->pmats,red->pmats,str);CHKERRQ(ierr);
1664b9ad928SBarry Smith   } else {
16790f1c854SHong Zhang     ierr = KSPSetOperators(red->ksp,pc->mat,pc->pmat,pc->flag);CHKERRQ(ierr);
1684b9ad928SBarry Smith   }
1690c24e6a1SHong Zhang   if (pc->setfromoptionscalled) {
1703e065800SHong Zhang     ierr = KSPSetFromOptions(red->ksp);CHKERRQ(ierr);
1710c24e6a1SHong Zhang   }
1723e065800SHong Zhang   ierr = KSPSetUp(red->ksp);CHKERRQ(ierr);
1734b9ad928SBarry Smith   PetscFunctionReturn(0);
1744b9ad928SBarry Smith }
1754b9ad928SBarry Smith 
1764b9ad928SBarry Smith #undef __FUNCT__
1774b9ad928SBarry Smith #define __FUNCT__ "PCApply_Redundant"
1786849ba73SBarry Smith static PetscErrorCode PCApply_Redundant(PC pc,Vec x,Vec y)
1794b9ad928SBarry Smith {
1804b9ad928SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
181dfbe8321SBarry Smith   PetscErrorCode ierr;
1823f457be1SHong Zhang   PetscScalar    *array;
1834b9ad928SBarry Smith 
1844b9ad928SBarry Smith   PetscFunctionBegin;
185*ddc54837SHong Zhang   if (!red->useparallelmat) {
186*ddc54837SHong Zhang     ierr = KSPSolve(red->ksp,x,y);CHKERRQ(ierr);
187*ddc54837SHong Zhang     PetscFunctionReturn(0);
188*ddc54837SHong Zhang   }
189*ddc54837SHong Zhang 
1903f457be1SHong Zhang   /* scatter x to xdup */
191ca9f406cSSatish Balay   ierr = VecScatterBegin(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
192ca9f406cSSatish Balay   ierr = VecScatterEnd(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1933f457be1SHong Zhang 
1943f457be1SHong Zhang   /* place xdup's local array into xsub */
1953f457be1SHong Zhang   ierr = VecGetArray(red->xdup,&array);CHKERRQ(ierr);
1963f457be1SHong Zhang   ierr = VecPlaceArray(red->xsub,(const PetscScalar*)array);CHKERRQ(ierr);
1974b9ad928SBarry Smith 
1984b9ad928SBarry Smith   /* apply preconditioner on each processor */
19983ab6a24SBarry Smith   ierr = KSPSolve(red->ksp,red->xsub,red->ysub);CHKERRQ(ierr);
2003f457be1SHong Zhang   ierr = VecResetArray(red->xsub);CHKERRQ(ierr);
2013f457be1SHong Zhang   ierr = VecRestoreArray(red->xdup,&array);CHKERRQ(ierr);
2024b9ad928SBarry Smith 
2033f457be1SHong Zhang   /* place ysub's local array into ydup */
2043f457be1SHong Zhang   ierr = VecGetArray(red->ysub,&array);CHKERRQ(ierr);
2053f457be1SHong Zhang   ierr = VecPlaceArray(red->ydup,(const PetscScalar*)array);CHKERRQ(ierr);
2063f457be1SHong Zhang 
2073f457be1SHong Zhang   /* scatter ydup to y */
208ca9f406cSSatish Balay   ierr = VecScatterBegin(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
209ca9f406cSSatish Balay   ierr = VecScatterEnd(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
2103f457be1SHong Zhang   ierr = VecResetArray(red->ydup);CHKERRQ(ierr);
2113f457be1SHong Zhang   ierr = VecRestoreArray(red->ysub,&array);CHKERRQ(ierr);
2124b9ad928SBarry Smith   PetscFunctionReturn(0);
2134b9ad928SBarry Smith }
2144b9ad928SBarry Smith 
2154b9ad928SBarry Smith #undef __FUNCT__
2161ea5a559SBarry Smith #define __FUNCT__ "PCReset_Redundant"
2171ea5a559SBarry Smith static PetscErrorCode PCReset_Redundant(PC pc)
2184b9ad928SBarry Smith {
2194b9ad928SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
220dfbe8321SBarry Smith   PetscErrorCode ierr;
2214b9ad928SBarry Smith 
2224b9ad928SBarry Smith   PetscFunctionBegin;
2236bf464f9SBarry Smith   ierr = VecScatterDestroy(&red->scatterin);CHKERRQ(ierr);
2246bf464f9SBarry Smith   ierr = VecScatterDestroy(&red->scatterout);CHKERRQ(ierr);
2256bf464f9SBarry Smith   ierr = VecDestroy(&red->ysub);CHKERRQ(ierr);
2266bf464f9SBarry Smith   ierr = VecDestroy(&red->xsub);CHKERRQ(ierr);
2276bf464f9SBarry Smith   ierr = VecDestroy(&red->xdup);CHKERRQ(ierr);
2286bf464f9SBarry Smith   ierr = VecDestroy(&red->ydup);CHKERRQ(ierr);
2296bf464f9SBarry Smith   ierr = MatDestroy(&red->pmats);CHKERRQ(ierr);
230546c30cbSLisandro Dalcin   if (red->ksp) {ierr = KSPReset(red->ksp);CHKERRQ(ierr);}
2311ea5a559SBarry Smith   PetscFunctionReturn(0);
2321ea5a559SBarry Smith }
2331ea5a559SBarry Smith 
2341ea5a559SBarry Smith #undef __FUNCT__
2351ea5a559SBarry Smith #define __FUNCT__ "PCDestroy_Redundant"
2361ea5a559SBarry Smith static PetscErrorCode PCDestroy_Redundant(PC pc)
2371ea5a559SBarry Smith {
2381ea5a559SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
2391ea5a559SBarry Smith   PetscErrorCode ierr;
2401ea5a559SBarry Smith 
2411ea5a559SBarry Smith   PetscFunctionBegin;
2421ea5a559SBarry Smith   ierr = PCReset_Redundant(pc);CHKERRQ(ierr);
2436bf464f9SBarry Smith   ierr = KSPDestroy(&red->ksp);CHKERRQ(ierr);
2446bf464f9SBarry Smith   ierr = PetscSubcommDestroy(&red->psubcomm);CHKERRQ(ierr);
245c31cb41cSBarry Smith   ierr = PetscFree(pc->data);CHKERRQ(ierr);
2464b9ad928SBarry Smith   PetscFunctionReturn(0);
2474b9ad928SBarry Smith }
2484b9ad928SBarry Smith 
2494b9ad928SBarry Smith #undef __FUNCT__
2504b9ad928SBarry Smith #define __FUNCT__ "PCSetFromOptions_Redundant"
2516849ba73SBarry Smith static PetscErrorCode PCSetFromOptions_Redundant(PC pc)
2524b9ad928SBarry Smith {
253a98ce0f4SHong Zhang   PetscErrorCode ierr;
254a98ce0f4SHong Zhang   PC_Redundant   *red = (PC_Redundant*)pc->data;
255a98ce0f4SHong Zhang 
2564b9ad928SBarry Smith   PetscFunctionBegin;
257a98ce0f4SHong Zhang   ierr = PetscOptionsHead("Redundant options");CHKERRQ(ierr);
25809a6bc64SHong Zhang   ierr = PetscOptionsInt("-pc_redundant_number","Number of redundant pc","PCRedundantSetNumber",red->nsubcomm,&red->nsubcomm,0);CHKERRQ(ierr);
259a98ce0f4SHong Zhang   ierr = PetscOptionsTail();CHKERRQ(ierr);
2604b9ad928SBarry Smith   PetscFunctionReturn(0);
2614b9ad928SBarry Smith }
2624b9ad928SBarry Smith 
2634b9ad928SBarry Smith #undef __FUNCT__
26409a6bc64SHong Zhang #define __FUNCT__ "PCRedundantSetNumber_Redundant"
265f7a08781SBarry Smith static PetscErrorCode  PCRedundantSetNumber_Redundant(PC pc,PetscInt nreds)
26609a6bc64SHong Zhang {
26709a6bc64SHong Zhang   PC_Redundant *red = (PC_Redundant*)pc->data;
26809a6bc64SHong Zhang 
26909a6bc64SHong Zhang   PetscFunctionBegin;
27009a6bc64SHong Zhang   red->nsubcomm = nreds;
27109a6bc64SHong Zhang   PetscFunctionReturn(0);
27209a6bc64SHong Zhang }
27309a6bc64SHong Zhang 
27409a6bc64SHong Zhang #undef __FUNCT__
27509a6bc64SHong Zhang #define __FUNCT__ "PCRedundantSetNumber"
27609a6bc64SHong Zhang /*@
27709a6bc64SHong Zhang    PCRedundantSetNumber - Sets the number of redundant preconditioner contexts.
27809a6bc64SHong Zhang 
2793f9fe445SBarry Smith    Logically Collective on PC
28009a6bc64SHong Zhang 
28109a6bc64SHong Zhang    Input Parameters:
28209a6bc64SHong Zhang +  pc - the preconditioner context
2839b21d695SBarry Smith -  nredundant - number of redundant preconditioner contexts; for example if you are using 64 MPI processes and
2849b21d695SBarry Smith                               use an nredundant of 4 there will be 4 parallel solves each on 16 = 64/4 processes.
28509a6bc64SHong Zhang 
28609a6bc64SHong Zhang    Level: advanced
28709a6bc64SHong Zhang 
28809a6bc64SHong Zhang .keywords: PC, redundant solve
28909a6bc64SHong Zhang @*/
2907087cfbeSBarry Smith PetscErrorCode  PCRedundantSetNumber(PC pc,PetscInt nredundant)
29109a6bc64SHong Zhang {
2924ac538c5SBarry Smith   PetscErrorCode ierr;
29309a6bc64SHong Zhang 
29409a6bc64SHong Zhang   PetscFunctionBegin;
2950700a824SBarry Smith   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
296ce94432eSBarry Smith   if (nredundant <= 0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "num of redundant pc %D must be positive",nredundant);
2974ac538c5SBarry Smith   ierr = PetscTryMethod(pc,"PCRedundantSetNumber_C",(PC,PetscInt),(pc,nredundant));CHKERRQ(ierr);
29809a6bc64SHong Zhang   PetscFunctionReturn(0);
29909a6bc64SHong Zhang }
30009a6bc64SHong Zhang 
30109a6bc64SHong Zhang #undef __FUNCT__
3024b9ad928SBarry Smith #define __FUNCT__ "PCRedundantSetScatter_Redundant"
303f7a08781SBarry Smith static PetscErrorCode  PCRedundantSetScatter_Redundant(PC pc,VecScatter in,VecScatter out)
3044b9ad928SBarry Smith {
3054b9ad928SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
306dfbe8321SBarry Smith   PetscErrorCode ierr;
3074b9ad928SBarry Smith 
3084b9ad928SBarry Smith   PetscFunctionBegin;
3094b9ad928SBarry Smith   ierr = PetscObjectReference((PetscObject)in);CHKERRQ(ierr);
3106bf464f9SBarry Smith   ierr = VecScatterDestroy(&red->scatterin);CHKERRQ(ierr);
3112fa5cd67SKarl Rupp 
312c3122656SLisandro Dalcin   red->scatterin  = in;
3132fa5cd67SKarl Rupp 
3144b9ad928SBarry Smith   ierr = PetscObjectReference((PetscObject)out);CHKERRQ(ierr);
3156bf464f9SBarry Smith   ierr = VecScatterDestroy(&red->scatterout);CHKERRQ(ierr);
316c3122656SLisandro Dalcin   red->scatterout = out;
3174b9ad928SBarry Smith   PetscFunctionReturn(0);
3184b9ad928SBarry Smith }
3194b9ad928SBarry Smith 
3204b9ad928SBarry Smith #undef __FUNCT__
3214b9ad928SBarry Smith #define __FUNCT__ "PCRedundantSetScatter"
3224b9ad928SBarry Smith /*@
3234b9ad928SBarry Smith    PCRedundantSetScatter - Sets the scatter used to copy values into the
3244b9ad928SBarry Smith      redundant local solve and the scatter to move them back into the global
3254b9ad928SBarry Smith      vector.
3264b9ad928SBarry Smith 
3273f9fe445SBarry Smith    Logically Collective on PC
3284b9ad928SBarry Smith 
3294b9ad928SBarry Smith    Input Parameters:
3304b9ad928SBarry Smith +  pc - the preconditioner context
3314b9ad928SBarry Smith .  in - the scatter to move the values in
3324b9ad928SBarry Smith -  out - the scatter to move them out
3334b9ad928SBarry Smith 
3344b9ad928SBarry Smith    Level: advanced
3354b9ad928SBarry Smith 
3364b9ad928SBarry Smith .keywords: PC, redundant solve
3374b9ad928SBarry Smith @*/
3387087cfbeSBarry Smith PetscErrorCode  PCRedundantSetScatter(PC pc,VecScatter in,VecScatter out)
3394b9ad928SBarry Smith {
3404ac538c5SBarry Smith   PetscErrorCode ierr;
3414b9ad928SBarry Smith 
3424b9ad928SBarry Smith   PetscFunctionBegin;
3430700a824SBarry Smith   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
3440700a824SBarry Smith   PetscValidHeaderSpecific(in,VEC_SCATTER_CLASSID,2);
3450700a824SBarry Smith   PetscValidHeaderSpecific(out,VEC_SCATTER_CLASSID,3);
3464ac538c5SBarry Smith   ierr = PetscTryMethod(pc,"PCRedundantSetScatter_C",(PC,VecScatter,VecScatter),(pc,in,out));CHKERRQ(ierr);
3474b9ad928SBarry Smith   PetscFunctionReturn(0);
3484b9ad928SBarry Smith }
3494b9ad928SBarry Smith 
3504b9ad928SBarry Smith #undef __FUNCT__
35183ab6a24SBarry Smith #define __FUNCT__ "PCRedundantGetKSP_Redundant"
352f7a08781SBarry Smith static PetscErrorCode  PCRedundantGetKSP_Redundant(PC pc,KSP *innerksp)
3534b9ad928SBarry Smith {
3545f06b7aaSBarry Smith   PetscErrorCode ierr;
3554b9ad928SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
3565f06b7aaSBarry Smith   MPI_Comm       comm,subcomm;
3575f06b7aaSBarry Smith   const char     *prefix;
3584b9ad928SBarry Smith 
3594b9ad928SBarry Smith   PetscFunctionBegin;
3605f06b7aaSBarry Smith   if (!red->psubcomm) {
3615f06b7aaSBarry Smith     ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
362d8a68f86SHong Zhang     ierr = PetscSubcommCreate(comm,&red->psubcomm);CHKERRQ(ierr);
363d8a68f86SHong Zhang     ierr = PetscSubcommSetNumber(red->psubcomm,red->nsubcomm);CHKERRQ(ierr);
364d8a68f86SHong Zhang     ierr = PetscSubcommSetType(red->psubcomm,PETSC_SUBCOMM_INTERLACED);CHKERRQ(ierr);
3655f06b7aaSBarry Smith     ierr = PetscLogObjectMemory(pc,sizeof(PetscSubcomm));CHKERRQ(ierr);
3665f06b7aaSBarry Smith 
3675f06b7aaSBarry Smith     /* create a new PC that processors in each subcomm have copy of */
3685f06b7aaSBarry Smith     subcomm = red->psubcomm->comm;
3692fa5cd67SKarl Rupp 
3705f06b7aaSBarry Smith     ierr = KSPCreate(subcomm,&red->ksp);CHKERRQ(ierr);
3715f06b7aaSBarry Smith     ierr = PetscObjectIncrementTabLevel((PetscObject)red->ksp,(PetscObject)pc,1);CHKERRQ(ierr);
3725f06b7aaSBarry Smith     ierr = PetscLogObjectParent(pc,red->ksp);CHKERRQ(ierr);
3735f06b7aaSBarry Smith     ierr = KSPSetType(red->ksp,KSPPREONLY);CHKERRQ(ierr);
3745f06b7aaSBarry Smith     ierr = KSPGetPC(red->ksp,&red->pc);CHKERRQ(ierr);
3755f06b7aaSBarry Smith     ierr = PCSetType(red->pc,PCLU);CHKERRQ(ierr);
3765f06b7aaSBarry Smith 
3775f06b7aaSBarry Smith     ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr);
3785f06b7aaSBarry Smith     ierr = KSPSetOptionsPrefix(red->ksp,prefix);CHKERRQ(ierr);
3795f06b7aaSBarry Smith     ierr = KSPAppendOptionsPrefix(red->ksp,"redundant_");CHKERRQ(ierr);
3805f06b7aaSBarry Smith   }
38183ab6a24SBarry Smith   *innerksp = red->ksp;
3824b9ad928SBarry Smith   PetscFunctionReturn(0);
3834b9ad928SBarry Smith }
3844b9ad928SBarry Smith 
3854b9ad928SBarry Smith #undef __FUNCT__
38683ab6a24SBarry Smith #define __FUNCT__ "PCRedundantGetKSP"
3874b9ad928SBarry Smith /*@
38883ab6a24SBarry Smith    PCRedundantGetKSP - Gets the less parallel KSP created by the redundant PC.
3894b9ad928SBarry Smith 
3904b9ad928SBarry Smith    Not Collective
3914b9ad928SBarry Smith 
3924b9ad928SBarry Smith    Input Parameter:
3934b9ad928SBarry Smith .  pc - the preconditioner context
3944b9ad928SBarry Smith 
3954b9ad928SBarry Smith    Output Parameter:
39683ab6a24SBarry Smith .  innerksp - the KSP on the smaller set of processes
3974b9ad928SBarry Smith 
3984b9ad928SBarry Smith    Level: advanced
3994b9ad928SBarry Smith 
4004b9ad928SBarry Smith .keywords: PC, redundant solve
4014b9ad928SBarry Smith @*/
40283ab6a24SBarry Smith PetscErrorCode  PCRedundantGetKSP(PC pc,KSP *innerksp)
4034b9ad928SBarry Smith {
4044ac538c5SBarry Smith   PetscErrorCode ierr;
4054b9ad928SBarry Smith 
4064b9ad928SBarry Smith   PetscFunctionBegin;
4070700a824SBarry Smith   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
40883ab6a24SBarry Smith   PetscValidPointer(innerksp,2);
40983ab6a24SBarry Smith   ierr = PetscTryMethod(pc,"PCRedundantGetKSP_C",(PC,KSP*),(pc,innerksp));CHKERRQ(ierr);
4104b9ad928SBarry Smith   PetscFunctionReturn(0);
4114b9ad928SBarry Smith }
4124b9ad928SBarry Smith 
4134b9ad928SBarry Smith #undef __FUNCT__
4144b9ad928SBarry Smith #define __FUNCT__ "PCRedundantGetOperators_Redundant"
415f7a08781SBarry Smith static PetscErrorCode  PCRedundantGetOperators_Redundant(PC pc,Mat *mat,Mat *pmat)
4164b9ad928SBarry Smith {
4174b9ad928SBarry Smith   PC_Redundant *red = (PC_Redundant*)pc->data;
4184b9ad928SBarry Smith 
4194b9ad928SBarry Smith   PetscFunctionBegin;
420b3804887SHong Zhang   if (mat)  *mat  = red->pmats;
421b3804887SHong Zhang   if (pmat) *pmat = red->pmats;
4224b9ad928SBarry Smith   PetscFunctionReturn(0);
4234b9ad928SBarry Smith }
4244b9ad928SBarry Smith 
4254b9ad928SBarry Smith #undef __FUNCT__
4264b9ad928SBarry Smith #define __FUNCT__ "PCRedundantGetOperators"
4274b9ad928SBarry Smith /*@
4284b9ad928SBarry Smith    PCRedundantGetOperators - gets the sequential matrix and preconditioner matrix
4294b9ad928SBarry Smith 
4304b9ad928SBarry Smith    Not Collective
4314b9ad928SBarry Smith 
4324b9ad928SBarry Smith    Input Parameter:
4334b9ad928SBarry Smith .  pc - the preconditioner context
4344b9ad928SBarry Smith 
4354b9ad928SBarry Smith    Output Parameters:
4364b9ad928SBarry Smith +  mat - the matrix
4374b9ad928SBarry Smith -  pmat - the (possibly different) preconditioner matrix
4384b9ad928SBarry Smith 
4394b9ad928SBarry Smith    Level: advanced
4404b9ad928SBarry Smith 
4414b9ad928SBarry Smith .keywords: PC, redundant solve
4424b9ad928SBarry Smith @*/
4437087cfbeSBarry Smith PetscErrorCode  PCRedundantGetOperators(PC pc,Mat *mat,Mat *pmat)
4444b9ad928SBarry Smith {
4454ac538c5SBarry Smith   PetscErrorCode ierr;
4464b9ad928SBarry Smith 
4474b9ad928SBarry Smith   PetscFunctionBegin;
4480700a824SBarry Smith   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
4494482741eSBarry Smith   if (mat)  PetscValidPointer(mat,2);
4504482741eSBarry Smith   if (pmat) PetscValidPointer(pmat,3);
4514ac538c5SBarry Smith   ierr = PetscTryMethod(pc,"PCRedundantGetOperators_C",(PC,Mat*,Mat*),(pc,mat,pmat));CHKERRQ(ierr);
4524b9ad928SBarry Smith   PetscFunctionReturn(0);
4534b9ad928SBarry Smith }
4544b9ad928SBarry Smith 
4554b9ad928SBarry Smith /* -------------------------------------------------------------------------------------*/
45637a17b4dSBarry Smith /*MC
45783ab6a24SBarry Smith      PCREDUNDANT - Runs a KSP solver with preconditioner for the entire problem on subgroups of processors
45837a17b4dSBarry Smith 
45983ab6a24SBarry Smith      Options for the redundant preconditioners can be set with -redundant_pc_xxx for the redundant KSP with -redundant_ksp_xxx
46037a17b4dSBarry Smith 
46109391456SBarry Smith   Options Database:
4629b21d695SBarry Smith .  -pc_redundant_number <n> - number of redundant solves, for example if you are using 64 MPI processes and
4639b21d695SBarry Smith                               use an n of 4 there will be 4 parallel solves each on 16 = 64/4 processes.
46409391456SBarry Smith 
46537a17b4dSBarry Smith    Level: intermediate
46637a17b4dSBarry Smith 
46783ab6a24SBarry Smith    Notes: The default KSP is preonly and the default PC is LU.
46883ab6a24SBarry Smith 
46983ab6a24SBarry Smith    Developer Notes: Note that PCSetInitialGuessNonzero()  is not used by this class but likely should be.
4709cfaa89bSBarry Smith 
47137a17b4dSBarry Smith .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PCRedundantSetScatter(),
47283ab6a24SBarry Smith            PCRedundantGetKSP(), PCRedundantGetOperators(), PCRedundantSetNumber()
47337a17b4dSBarry Smith M*/
47437a17b4dSBarry Smith 
4754b9ad928SBarry Smith #undef __FUNCT__
4764b9ad928SBarry Smith #define __FUNCT__ "PCCreate_Redundant"
4778cc058d9SJed Brown PETSC_EXTERN PetscErrorCode PCCreate_Redundant(PC pc)
4784b9ad928SBarry Smith {
479dfbe8321SBarry Smith   PetscErrorCode ierr;
4804b9ad928SBarry Smith   PC_Redundant   *red;
48169db28dcSHong Zhang   PetscMPIInt    size;
4823f457be1SHong Zhang 
4834b9ad928SBarry Smith   PetscFunctionBegin;
48438f2d2fdSLisandro Dalcin   ierr = PetscNewLog(pc,PC_Redundant,&red);CHKERRQ(ierr);
485ce94432eSBarry Smith   ierr = MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);CHKERRQ(ierr);
4862fa5cd67SKarl Rupp 
48769db28dcSHong Zhang   red->nsubcomm       = size;
4884b9ad928SBarry Smith   red->useparallelmat = PETSC_TRUE;
4891fbd8f88SHong Zhang   pc->data            = (void*)red;
4904b9ad928SBarry Smith 
4914b9ad928SBarry Smith   pc->ops->apply          = PCApply_Redundant;
4924b9ad928SBarry Smith   pc->ops->applytranspose = 0;
4934b9ad928SBarry Smith   pc->ops->setup          = PCSetUp_Redundant;
4944b9ad928SBarry Smith   pc->ops->destroy        = PCDestroy_Redundant;
4951ea5a559SBarry Smith   pc->ops->reset          = PCReset_Redundant;
4964b9ad928SBarry Smith   pc->ops->setfromoptions = PCSetFromOptions_Redundant;
4974b9ad928SBarry Smith   pc->ops->view           = PCView_Redundant;
4982fa5cd67SKarl Rupp 
499bdf89e91SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCRedundantSetScatter_C",PCRedundantSetScatter_Redundant);CHKERRQ(ierr);
500bdf89e91SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCRedundantSetNumber_C",PCRedundantSetNumber_Redundant);CHKERRQ(ierr);
501bdf89e91SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCRedundantGetKSP_C",PCRedundantGetKSP_Redundant);CHKERRQ(ierr);
502bdf89e91SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCRedundantGetOperators_C",PCRedundantGetOperators_Redundant);CHKERRQ(ierr);
5034b9ad928SBarry Smith   PetscFunctionReturn(0);
5044b9ad928SBarry Smith }
505b2573a8aSBarry Smith 
506