xref: /petsc/src/ksp/pc/impls/redundant/redundant.c (revision f68be91c59f5f8ea02c6bea4ec423ba983a84562)
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 
51e37c6257SHong Zhang extern PetscErrorCode MatGetRedundantMatrix_MPIAIJ_psubcomm(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;
591b81debcSHong Zhang   PetscInt       mstart,mend,mlocal,M;
6013f74950SBarry Smith   PetscMPIInt    size;
61ce94432eSBarry Smith   MPI_Comm       comm,subcomm;
62ddc54837SHong Zhang   Vec            x;
631fbd8f88SHong Zhang   const char     *prefix;
643f457be1SHong Zhang 
654b9ad928SBarry Smith   PetscFunctionBegin;
66ce94432eSBarry Smith   ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
67ddc54837SHong Zhang 
68ddc54837SHong Zhang   /* if pmatrix set by user is sequential then we do not need to gather the parallel matrix */
69ddc54837SHong Zhang   ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
70ddc54837SHong Zhang   if (size == 1) red->useparallelmat = PETSC_FALSE;
711fbd8f88SHong Zhang 
724b9ad928SBarry Smith   if (!pc->setupcalled) {
731b81debcSHong Zhang     PetscInt mloc_sub;
745f06b7aaSBarry Smith     if (!red->psubcomm) {
75d8a68f86SHong Zhang       ierr = PetscSubcommCreate(comm,&red->psubcomm);CHKERRQ(ierr);
76d8a68f86SHong Zhang       ierr = PetscSubcommSetNumber(red->psubcomm,red->nsubcomm);CHKERRQ(ierr);
77d8a68f86SHong Zhang       ierr = PetscSubcommSetType(red->psubcomm,PETSC_SUBCOMM_INTERLACED);CHKERRQ(ierr);
78*f68be91cSHong Zhang       /* enable runtime switch of psubcomm type, e.g., '-psubcomm_type contiguous */
79*f68be91cSHong Zhang       ierr = PetscSubcommSetFromOptions(red->psubcomm);CHKERRQ(ierr);
80*f68be91cSHong Zhang 
811fbd8f88SHong Zhang       ierr = PetscLogObjectMemory(pc,sizeof(PetscSubcomm));CHKERRQ(ierr);
821fbd8f88SHong Zhang 
831fbd8f88SHong Zhang       /* create a new PC that processors in each subcomm have copy of */
840d7810c8SBarry Smith       subcomm = red->psubcomm->comm;
852fa5cd67SKarl Rupp 
865f06b7aaSBarry Smith       ierr = KSPCreate(subcomm,&red->ksp);CHKERRQ(ierr);
875f06b7aaSBarry Smith       ierr = PetscObjectIncrementTabLevel((PetscObject)red->ksp,(PetscObject)pc,1);CHKERRQ(ierr);
885f06b7aaSBarry Smith       ierr = PetscLogObjectParent(pc,red->ksp);CHKERRQ(ierr);
895f06b7aaSBarry Smith       ierr = KSPSetType(red->ksp,KSPPREONLY);CHKERRQ(ierr);
905f06b7aaSBarry Smith       ierr = KSPGetPC(red->ksp,&red->pc);CHKERRQ(ierr);
91cf52b8b1SHong Zhang       ierr = PCSetType(red->pc,PCLU);CHKERRQ(ierr);
92cf52b8b1SHong Zhang 
931fbd8f88SHong Zhang       ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr);
945f06b7aaSBarry Smith       ierr = KSPSetOptionsPrefix(red->ksp,prefix);CHKERRQ(ierr);
955f06b7aaSBarry Smith       ierr = KSPAppendOptionsPrefix(red->ksp,"redundant_");CHKERRQ(ierr);
961b81debcSHong Zhang     } else {
971b81debcSHong Zhang       subcomm = red->psubcomm->comm;
981b81debcSHong Zhang     }
991fbd8f88SHong Zhang 
1001b81debcSHong Zhang     if (red->useparallelmat) {
1011b81debcSHong Zhang       /* grab the parallel matrix and put it into processors of a subcomminicator */
102e37c6257SHong Zhang       ierr = MatGetRedundantMatrix_MPIAIJ_psubcomm(pc->pmat,red->psubcomm->n,red->psubcomm,MAT_INITIAL_MATRIX,&red->pmats);CHKERRQ(ierr);
1031b81debcSHong Zhang       ierr = KSPSetOperators(red->ksp,red->pmats,red->pmats,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
1044b9ad928SBarry Smith 
1051b81debcSHong Zhang       /* get working vectors xsub and ysub */
1061b81debcSHong Zhang       ierr = MatGetVecs(red->pmats,&red->xsub,&red->ysub);CHKERRQ(ierr);
1072fa5cd67SKarl Rupp 
1088b96b0d2SHong Zhang       /* create working vectors xdup and ydup.
1098b96b0d2SHong Zhang        xdup concatenates all xsub's contigously to form a mpi vector over dupcomm  (see PetscSubcommCreate_interlaced())
1108b96b0d2SHong Zhang        ydup concatenates all ysub and has empty local arrays because ysub's arrays will be place into it.
111ce94432eSBarry Smith        Note: we use communicator dupcomm, not PetscObjectComm((PetscObject)pc)! */
1121b81debcSHong Zhang       ierr = MatGetLocalSize(red->pmats,&mloc_sub,NULL);CHKERRQ(ierr);
1131fbd8f88SHong Zhang       ierr = VecCreateMPI(red->psubcomm->dupparent,mloc_sub,PETSC_DECIDE,&red->xdup);CHKERRQ(ierr);
1140298fd71SBarry Smith       ierr = VecCreateMPIWithArray(red->psubcomm->dupparent,1,mloc_sub,PETSC_DECIDE,NULL,&red->ydup);CHKERRQ(ierr);
1153f457be1SHong Zhang 
116*f68be91cSHong Zhang       /* create vecscatters */
117*f68be91cSHong Zhang       if (!red->scatterin) { /* efficiency of scatterin is independent from psubcomm_type! */
1183f457be1SHong Zhang         IS       is1,is2;
1193f457be1SHong Zhang         PetscInt *idx1,*idx2,i,j,k;
12045fc02eaSBarry Smith 
1211b81debcSHong Zhang         ierr = MatGetVecs(pc->pmat,&x,0);CHKERRQ(ierr);
1221b81debcSHong Zhang         ierr = VecGetSize(x,&M);CHKERRQ(ierr);
1231b81debcSHong Zhang         ierr = VecGetOwnershipRange(x,&mstart,&mend);CHKERRQ(ierr);
1241b81debcSHong Zhang         mlocal = mend - mstart;
1251714dc9eSHong Zhang         ierr = PetscMalloc2(red->psubcomm->n*mlocal,PetscInt,&idx1,red->psubcomm->n*mlocal,PetscInt,&idx2);CHKERRQ(ierr);
1263f457be1SHong Zhang         j    = 0;
1271fbd8f88SHong Zhang         for (k=0; k<red->psubcomm->n; k++) {
1283f457be1SHong Zhang           for (i=mstart; i<mend; i++) {
1293f457be1SHong Zhang             idx1[j]   = i;
130ddc54837SHong Zhang             idx2[j++] = i + M*k;
1313f457be1SHong Zhang           }
1323f457be1SHong Zhang         }
13370b3c8c7SBarry Smith         ierr = ISCreateGeneral(comm,red->psubcomm->n*mlocal,idx1,PETSC_COPY_VALUES,&is1);CHKERRQ(ierr);
13470b3c8c7SBarry Smith         ierr = ISCreateGeneral(comm,red->psubcomm->n*mlocal,idx2,PETSC_COPY_VALUES,&is2);CHKERRQ(ierr);
135ddc54837SHong Zhang         ierr = VecScatterCreate(x,is1,red->xdup,is2,&red->scatterin);CHKERRQ(ierr);
136fcfd50ebSBarry Smith         ierr = ISDestroy(&is1);CHKERRQ(ierr);
137fcfd50ebSBarry Smith         ierr = ISDestroy(&is2);CHKERRQ(ierr);
1383f457be1SHong Zhang 
139*f68be91cSHong Zhang         /* efficiency of scatterout depends on psubcomm_type! Impl below is good for PETSC_SUBCOMM_INTERLACED */
140ddc54837SHong Zhang         ierr = ISCreateStride(comm,mlocal,mstart+ red->psubcomm->color*M,1,&is1);CHKERRQ(ierr);
1413f457be1SHong Zhang         ierr = ISCreateStride(comm,mlocal,mstart,1,&is2);CHKERRQ(ierr);
142ddc54837SHong Zhang         ierr = VecScatterCreate(red->xdup,is1,x,is2,&red->scatterout);CHKERRQ(ierr);
143fcfd50ebSBarry Smith         ierr = ISDestroy(&is1);CHKERRQ(ierr);
144fcfd50ebSBarry Smith         ierr = ISDestroy(&is2);CHKERRQ(ierr);
1451d79065fSBarry Smith         ierr = PetscFree2(idx1,idx2);CHKERRQ(ierr);
146ddc54837SHong Zhang         ierr = VecDestroy(&x);CHKERRQ(ierr);
1471b81debcSHong Zhang       }
148ab661555SHong Zhang     } else { /* !red->useparallelmat */
149ab661555SHong Zhang       ierr = KSPSetOperators(red->ksp,pc->mat,pc->pmat,pc->flag);CHKERRQ(ierr);
1501b81debcSHong Zhang     }
151ab661555SHong Zhang   } else { /* pc->setupcalled */
1524b9ad928SBarry Smith     if (red->useparallelmat) {
153ab661555SHong Zhang       MatReuse       reuse;
1541b81debcSHong Zhang       /* grab the parallel matrix and put it into processors of a subcomminicator */
1551b81debcSHong Zhang       /*--------------------------------------------------------------------------*/
156ab661555SHong Zhang       if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1574b9ad928SBarry Smith         /* destroy old matrices */
1586bf464f9SBarry Smith         ierr  = MatDestroy(&red->pmats);CHKERRQ(ierr);
159ab661555SHong Zhang         reuse = MAT_INITIAL_MATRIX;
1604b9ad928SBarry Smith       } else {
161ab661555SHong Zhang         reuse = MAT_REUSE_MATRIX;
162ab661555SHong Zhang       }
163e37c6257SHong Zhang       ierr = MatGetRedundantMatrix_MPIAIJ_psubcomm(pc->pmat,red->psubcomm->n,red->psubcomm,reuse,&red->pmats);CHKERRQ(ierr);
164ab661555SHong Zhang       ierr = KSPSetOperators(red->ksp,red->pmats,red->pmats,pc->flag);CHKERRQ(ierr);
165ab661555SHong Zhang     } else { /* !red->useparallelmat */
16690f1c854SHong Zhang       ierr = KSPSetOperators(red->ksp,pc->mat,pc->pmat,pc->flag);CHKERRQ(ierr);
1674b9ad928SBarry Smith     }
168ab661555SHong Zhang   }
1691b81debcSHong Zhang 
1700c24e6a1SHong Zhang   if (pc->setfromoptionscalled) {
1713e065800SHong Zhang     ierr = KSPSetFromOptions(red->ksp);CHKERRQ(ierr);
1720c24e6a1SHong Zhang   }
1733e065800SHong Zhang   ierr = KSPSetUp(red->ksp);CHKERRQ(ierr);
1744b9ad928SBarry Smith   PetscFunctionReturn(0);
1754b9ad928SBarry Smith }
1764b9ad928SBarry Smith 
1774b9ad928SBarry Smith #undef __FUNCT__
1784b9ad928SBarry Smith #define __FUNCT__ "PCApply_Redundant"
1796849ba73SBarry Smith static PetscErrorCode PCApply_Redundant(PC pc,Vec x,Vec y)
1804b9ad928SBarry Smith {
1814b9ad928SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
182dfbe8321SBarry Smith   PetscErrorCode ierr;
1833f457be1SHong Zhang   PetscScalar    *array;
1844b9ad928SBarry Smith 
1854b9ad928SBarry Smith   PetscFunctionBegin;
186ddc54837SHong Zhang   if (!red->useparallelmat) {
187ddc54837SHong Zhang     ierr = KSPSolve(red->ksp,x,y);CHKERRQ(ierr);
188ddc54837SHong Zhang     PetscFunctionReturn(0);
189ddc54837SHong Zhang   }
190ddc54837SHong Zhang 
1913f457be1SHong Zhang   /* scatter x to xdup */
192ca9f406cSSatish Balay   ierr = VecScatterBegin(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
193ca9f406cSSatish Balay   ierr = VecScatterEnd(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1943f457be1SHong Zhang 
1953f457be1SHong Zhang   /* place xdup's local array into xsub */
1963f457be1SHong Zhang   ierr = VecGetArray(red->xdup,&array);CHKERRQ(ierr);
1973f457be1SHong Zhang   ierr = VecPlaceArray(red->xsub,(const PetscScalar*)array);CHKERRQ(ierr);
1984b9ad928SBarry Smith 
1994b9ad928SBarry Smith   /* apply preconditioner on each processor */
20083ab6a24SBarry Smith   ierr = KSPSolve(red->ksp,red->xsub,red->ysub);CHKERRQ(ierr);
2013f457be1SHong Zhang   ierr = VecResetArray(red->xsub);CHKERRQ(ierr);
2023f457be1SHong Zhang   ierr = VecRestoreArray(red->xdup,&array);CHKERRQ(ierr);
2034b9ad928SBarry Smith 
2043f457be1SHong Zhang   /* place ysub's local array into ydup */
2053f457be1SHong Zhang   ierr = VecGetArray(red->ysub,&array);CHKERRQ(ierr);
2063f457be1SHong Zhang   ierr = VecPlaceArray(red->ydup,(const PetscScalar*)array);CHKERRQ(ierr);
2073f457be1SHong Zhang 
2083f457be1SHong Zhang   /* scatter ydup to y */
209ca9f406cSSatish Balay   ierr = VecScatterBegin(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
210ca9f406cSSatish Balay   ierr = VecScatterEnd(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
2113f457be1SHong Zhang   ierr = VecResetArray(red->ydup);CHKERRQ(ierr);
2123f457be1SHong Zhang   ierr = VecRestoreArray(red->ysub,&array);CHKERRQ(ierr);
2134b9ad928SBarry Smith   PetscFunctionReturn(0);
2144b9ad928SBarry Smith }
2154b9ad928SBarry Smith 
2164b9ad928SBarry Smith #undef __FUNCT__
2171ea5a559SBarry Smith #define __FUNCT__ "PCReset_Redundant"
2181ea5a559SBarry Smith static PetscErrorCode PCReset_Redundant(PC pc)
2194b9ad928SBarry Smith {
2204b9ad928SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
221dfbe8321SBarry Smith   PetscErrorCode ierr;
2224b9ad928SBarry Smith 
2234b9ad928SBarry Smith   PetscFunctionBegin;
2241b81debcSHong Zhang   if (red->useparallelmat) {
2256bf464f9SBarry Smith     ierr = VecScatterDestroy(&red->scatterin);CHKERRQ(ierr);
2266bf464f9SBarry Smith     ierr = VecScatterDestroy(&red->scatterout);CHKERRQ(ierr);
2276bf464f9SBarry Smith     ierr = VecDestroy(&red->ysub);CHKERRQ(ierr);
2286bf464f9SBarry Smith     ierr = VecDestroy(&red->xsub);CHKERRQ(ierr);
2296bf464f9SBarry Smith     ierr = VecDestroy(&red->xdup);CHKERRQ(ierr);
2306bf464f9SBarry Smith     ierr = VecDestroy(&red->ydup);CHKERRQ(ierr);
2311b81debcSHong Zhang   }
2326bf464f9SBarry Smith   ierr = MatDestroy(&red->pmats);CHKERRQ(ierr);
2331b81debcSHong Zhang   ierr = KSPReset(red->ksp);CHKERRQ(ierr);
2341ea5a559SBarry Smith   PetscFunctionReturn(0);
2351ea5a559SBarry Smith }
2361ea5a559SBarry Smith 
2371ea5a559SBarry Smith #undef __FUNCT__
2381ea5a559SBarry Smith #define __FUNCT__ "PCDestroy_Redundant"
2391ea5a559SBarry Smith static PetscErrorCode PCDestroy_Redundant(PC pc)
2401ea5a559SBarry Smith {
2411ea5a559SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
2421ea5a559SBarry Smith   PetscErrorCode ierr;
2431ea5a559SBarry Smith 
2441ea5a559SBarry Smith   PetscFunctionBegin;
2451ea5a559SBarry Smith   ierr = PCReset_Redundant(pc);CHKERRQ(ierr);
2466bf464f9SBarry Smith   ierr = KSPDestroy(&red->ksp);CHKERRQ(ierr);
2476bf464f9SBarry Smith   ierr = PetscSubcommDestroy(&red->psubcomm);CHKERRQ(ierr);
248c31cb41cSBarry Smith   ierr = PetscFree(pc->data);CHKERRQ(ierr);
2494b9ad928SBarry Smith   PetscFunctionReturn(0);
2504b9ad928SBarry Smith }
2514b9ad928SBarry Smith 
2524b9ad928SBarry Smith #undef __FUNCT__
2534b9ad928SBarry Smith #define __FUNCT__ "PCSetFromOptions_Redundant"
2546849ba73SBarry Smith static PetscErrorCode PCSetFromOptions_Redundant(PC pc)
2554b9ad928SBarry Smith {
256a98ce0f4SHong Zhang   PetscErrorCode ierr;
257a98ce0f4SHong Zhang   PC_Redundant   *red = (PC_Redundant*)pc->data;
258a98ce0f4SHong Zhang 
2594b9ad928SBarry Smith   PetscFunctionBegin;
260a98ce0f4SHong Zhang   ierr = PetscOptionsHead("Redundant options");CHKERRQ(ierr);
26109a6bc64SHong Zhang   ierr = PetscOptionsInt("-pc_redundant_number","Number of redundant pc","PCRedundantSetNumber",red->nsubcomm,&red->nsubcomm,0);CHKERRQ(ierr);
262a98ce0f4SHong Zhang   ierr = PetscOptionsTail();CHKERRQ(ierr);
2634b9ad928SBarry Smith   PetscFunctionReturn(0);
2644b9ad928SBarry Smith }
2654b9ad928SBarry Smith 
2664b9ad928SBarry Smith #undef __FUNCT__
26709a6bc64SHong Zhang #define __FUNCT__ "PCRedundantSetNumber_Redundant"
268f7a08781SBarry Smith static PetscErrorCode  PCRedundantSetNumber_Redundant(PC pc,PetscInt nreds)
26909a6bc64SHong Zhang {
27009a6bc64SHong Zhang   PC_Redundant *red = (PC_Redundant*)pc->data;
27109a6bc64SHong Zhang 
27209a6bc64SHong Zhang   PetscFunctionBegin;
27309a6bc64SHong Zhang   red->nsubcomm = nreds;
27409a6bc64SHong Zhang   PetscFunctionReturn(0);
27509a6bc64SHong Zhang }
27609a6bc64SHong Zhang 
27709a6bc64SHong Zhang #undef __FUNCT__
27809a6bc64SHong Zhang #define __FUNCT__ "PCRedundantSetNumber"
27909a6bc64SHong Zhang /*@
28009a6bc64SHong Zhang    PCRedundantSetNumber - Sets the number of redundant preconditioner contexts.
28109a6bc64SHong Zhang 
2823f9fe445SBarry Smith    Logically Collective on PC
28309a6bc64SHong Zhang 
28409a6bc64SHong Zhang    Input Parameters:
28509a6bc64SHong Zhang +  pc - the preconditioner context
2869b21d695SBarry Smith -  nredundant - number of redundant preconditioner contexts; for example if you are using 64 MPI processes and
2879b21d695SBarry Smith                               use an nredundant of 4 there will be 4 parallel solves each on 16 = 64/4 processes.
28809a6bc64SHong Zhang 
28909a6bc64SHong Zhang    Level: advanced
29009a6bc64SHong Zhang 
29109a6bc64SHong Zhang .keywords: PC, redundant solve
29209a6bc64SHong Zhang @*/
2937087cfbeSBarry Smith PetscErrorCode  PCRedundantSetNumber(PC pc,PetscInt nredundant)
29409a6bc64SHong Zhang {
2954ac538c5SBarry Smith   PetscErrorCode ierr;
29609a6bc64SHong Zhang 
29709a6bc64SHong Zhang   PetscFunctionBegin;
2980700a824SBarry Smith   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
299ce94432eSBarry Smith   if (nredundant <= 0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "num of redundant pc %D must be positive",nredundant);
3004ac538c5SBarry Smith   ierr = PetscTryMethod(pc,"PCRedundantSetNumber_C",(PC,PetscInt),(pc,nredundant));CHKERRQ(ierr);
30109a6bc64SHong Zhang   PetscFunctionReturn(0);
30209a6bc64SHong Zhang }
30309a6bc64SHong Zhang 
30409a6bc64SHong Zhang #undef __FUNCT__
3054b9ad928SBarry Smith #define __FUNCT__ "PCRedundantSetScatter_Redundant"
306f7a08781SBarry Smith static PetscErrorCode  PCRedundantSetScatter_Redundant(PC pc,VecScatter in,VecScatter out)
3074b9ad928SBarry Smith {
3084b9ad928SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
309dfbe8321SBarry Smith   PetscErrorCode ierr;
3104b9ad928SBarry Smith 
3114b9ad928SBarry Smith   PetscFunctionBegin;
3124b9ad928SBarry Smith   ierr = PetscObjectReference((PetscObject)in);CHKERRQ(ierr);
3136bf464f9SBarry Smith   ierr = VecScatterDestroy(&red->scatterin);CHKERRQ(ierr);
3142fa5cd67SKarl Rupp 
315c3122656SLisandro Dalcin   red->scatterin  = in;
3162fa5cd67SKarl Rupp 
3174b9ad928SBarry Smith   ierr = PetscObjectReference((PetscObject)out);CHKERRQ(ierr);
3186bf464f9SBarry Smith   ierr = VecScatterDestroy(&red->scatterout);CHKERRQ(ierr);
319c3122656SLisandro Dalcin   red->scatterout = out;
3204b9ad928SBarry Smith   PetscFunctionReturn(0);
3214b9ad928SBarry Smith }
3224b9ad928SBarry Smith 
3234b9ad928SBarry Smith #undef __FUNCT__
3244b9ad928SBarry Smith #define __FUNCT__ "PCRedundantSetScatter"
3254b9ad928SBarry Smith /*@
3264b9ad928SBarry Smith    PCRedundantSetScatter - Sets the scatter used to copy values into the
3274b9ad928SBarry Smith      redundant local solve and the scatter to move them back into the global
3284b9ad928SBarry Smith      vector.
3294b9ad928SBarry Smith 
3303f9fe445SBarry Smith    Logically Collective on PC
3314b9ad928SBarry Smith 
3324b9ad928SBarry Smith    Input Parameters:
3334b9ad928SBarry Smith +  pc - the preconditioner context
3344b9ad928SBarry Smith .  in - the scatter to move the values in
3354b9ad928SBarry Smith -  out - the scatter to move them out
3364b9ad928SBarry Smith 
3374b9ad928SBarry Smith    Level: advanced
3384b9ad928SBarry Smith 
3394b9ad928SBarry Smith .keywords: PC, redundant solve
3404b9ad928SBarry Smith @*/
3417087cfbeSBarry Smith PetscErrorCode  PCRedundantSetScatter(PC pc,VecScatter in,VecScatter out)
3424b9ad928SBarry Smith {
3434ac538c5SBarry Smith   PetscErrorCode ierr;
3444b9ad928SBarry Smith 
3454b9ad928SBarry Smith   PetscFunctionBegin;
3460700a824SBarry Smith   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
3470700a824SBarry Smith   PetscValidHeaderSpecific(in,VEC_SCATTER_CLASSID,2);
3480700a824SBarry Smith   PetscValidHeaderSpecific(out,VEC_SCATTER_CLASSID,3);
3494ac538c5SBarry Smith   ierr = PetscTryMethod(pc,"PCRedundantSetScatter_C",(PC,VecScatter,VecScatter),(pc,in,out));CHKERRQ(ierr);
3504b9ad928SBarry Smith   PetscFunctionReturn(0);
3514b9ad928SBarry Smith }
3524b9ad928SBarry Smith 
3534b9ad928SBarry Smith #undef __FUNCT__
35483ab6a24SBarry Smith #define __FUNCT__ "PCRedundantGetKSP_Redundant"
355f7a08781SBarry Smith static PetscErrorCode  PCRedundantGetKSP_Redundant(PC pc,KSP *innerksp)
3564b9ad928SBarry Smith {
3575f06b7aaSBarry Smith   PetscErrorCode ierr;
3584b9ad928SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
3595f06b7aaSBarry Smith   MPI_Comm       comm,subcomm;
3605f06b7aaSBarry Smith   const char     *prefix;
3614b9ad928SBarry Smith 
3624b9ad928SBarry Smith   PetscFunctionBegin;
3635f06b7aaSBarry Smith   if (!red->psubcomm) {
3645f06b7aaSBarry Smith     ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
365d8a68f86SHong Zhang     ierr = PetscSubcommCreate(comm,&red->psubcomm);CHKERRQ(ierr);
366d8a68f86SHong Zhang     ierr = PetscSubcommSetNumber(red->psubcomm,red->nsubcomm);CHKERRQ(ierr);
367d8a68f86SHong Zhang     ierr = PetscSubcommSetType(red->psubcomm,PETSC_SUBCOMM_INTERLACED);CHKERRQ(ierr);
3685f06b7aaSBarry Smith     ierr = PetscLogObjectMemory(pc,sizeof(PetscSubcomm));CHKERRQ(ierr);
3695f06b7aaSBarry Smith 
3705f06b7aaSBarry Smith     /* create a new PC that processors in each subcomm have copy of */
3715f06b7aaSBarry Smith     subcomm = red->psubcomm->comm;
3722fa5cd67SKarl Rupp 
3735f06b7aaSBarry Smith     ierr = KSPCreate(subcomm,&red->ksp);CHKERRQ(ierr);
3745f06b7aaSBarry Smith     ierr = PetscObjectIncrementTabLevel((PetscObject)red->ksp,(PetscObject)pc,1);CHKERRQ(ierr);
3755f06b7aaSBarry Smith     ierr = PetscLogObjectParent(pc,red->ksp);CHKERRQ(ierr);
3765f06b7aaSBarry Smith     ierr = KSPSetType(red->ksp,KSPPREONLY);CHKERRQ(ierr);
3775f06b7aaSBarry Smith     ierr = KSPGetPC(red->ksp,&red->pc);CHKERRQ(ierr);
3785f06b7aaSBarry Smith     ierr = PCSetType(red->pc,PCLU);CHKERRQ(ierr);
3795f06b7aaSBarry Smith 
3805f06b7aaSBarry Smith     ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr);
3815f06b7aaSBarry Smith     ierr = KSPSetOptionsPrefix(red->ksp,prefix);CHKERRQ(ierr);
3825f06b7aaSBarry Smith     ierr = KSPAppendOptionsPrefix(red->ksp,"redundant_");CHKERRQ(ierr);
3835f06b7aaSBarry Smith   }
38483ab6a24SBarry Smith   *innerksp = red->ksp;
3854b9ad928SBarry Smith   PetscFunctionReturn(0);
3864b9ad928SBarry Smith }
3874b9ad928SBarry Smith 
3884b9ad928SBarry Smith #undef __FUNCT__
38983ab6a24SBarry Smith #define __FUNCT__ "PCRedundantGetKSP"
3904b9ad928SBarry Smith /*@
39183ab6a24SBarry Smith    PCRedundantGetKSP - Gets the less parallel KSP created by the redundant PC.
3924b9ad928SBarry Smith 
3934b9ad928SBarry Smith    Not Collective
3944b9ad928SBarry Smith 
3954b9ad928SBarry Smith    Input Parameter:
3964b9ad928SBarry Smith .  pc - the preconditioner context
3974b9ad928SBarry Smith 
3984b9ad928SBarry Smith    Output Parameter:
39983ab6a24SBarry Smith .  innerksp - the KSP on the smaller set of processes
4004b9ad928SBarry Smith 
4014b9ad928SBarry Smith    Level: advanced
4024b9ad928SBarry Smith 
4034b9ad928SBarry Smith .keywords: PC, redundant solve
4044b9ad928SBarry Smith @*/
40583ab6a24SBarry Smith PetscErrorCode  PCRedundantGetKSP(PC pc,KSP *innerksp)
4064b9ad928SBarry Smith {
4074ac538c5SBarry Smith   PetscErrorCode ierr;
4084b9ad928SBarry Smith 
4094b9ad928SBarry Smith   PetscFunctionBegin;
4100700a824SBarry Smith   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
41183ab6a24SBarry Smith   PetscValidPointer(innerksp,2);
41283ab6a24SBarry Smith   ierr = PetscTryMethod(pc,"PCRedundantGetKSP_C",(PC,KSP*),(pc,innerksp));CHKERRQ(ierr);
4134b9ad928SBarry Smith   PetscFunctionReturn(0);
4144b9ad928SBarry Smith }
4154b9ad928SBarry Smith 
4164b9ad928SBarry Smith #undef __FUNCT__
4174b9ad928SBarry Smith #define __FUNCT__ "PCRedundantGetOperators_Redundant"
418f7a08781SBarry Smith static PetscErrorCode  PCRedundantGetOperators_Redundant(PC pc,Mat *mat,Mat *pmat)
4194b9ad928SBarry Smith {
4204b9ad928SBarry Smith   PC_Redundant *red = (PC_Redundant*)pc->data;
4214b9ad928SBarry Smith 
4224b9ad928SBarry Smith   PetscFunctionBegin;
423b3804887SHong Zhang   if (mat)  *mat  = red->pmats;
424b3804887SHong Zhang   if (pmat) *pmat = red->pmats;
4254b9ad928SBarry Smith   PetscFunctionReturn(0);
4264b9ad928SBarry Smith }
4274b9ad928SBarry Smith 
4284b9ad928SBarry Smith #undef __FUNCT__
4294b9ad928SBarry Smith #define __FUNCT__ "PCRedundantGetOperators"
4304b9ad928SBarry Smith /*@
4314b9ad928SBarry Smith    PCRedundantGetOperators - gets the sequential matrix and preconditioner matrix
4324b9ad928SBarry Smith 
4334b9ad928SBarry Smith    Not Collective
4344b9ad928SBarry Smith 
4354b9ad928SBarry Smith    Input Parameter:
4364b9ad928SBarry Smith .  pc - the preconditioner context
4374b9ad928SBarry Smith 
4384b9ad928SBarry Smith    Output Parameters:
4394b9ad928SBarry Smith +  mat - the matrix
4404b9ad928SBarry Smith -  pmat - the (possibly different) preconditioner matrix
4414b9ad928SBarry Smith 
4424b9ad928SBarry Smith    Level: advanced
4434b9ad928SBarry Smith 
4444b9ad928SBarry Smith .keywords: PC, redundant solve
4454b9ad928SBarry Smith @*/
4467087cfbeSBarry Smith PetscErrorCode  PCRedundantGetOperators(PC pc,Mat *mat,Mat *pmat)
4474b9ad928SBarry Smith {
4484ac538c5SBarry Smith   PetscErrorCode ierr;
4494b9ad928SBarry Smith 
4504b9ad928SBarry Smith   PetscFunctionBegin;
4510700a824SBarry Smith   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
4524482741eSBarry Smith   if (mat)  PetscValidPointer(mat,2);
4534482741eSBarry Smith   if (pmat) PetscValidPointer(pmat,3);
4544ac538c5SBarry Smith   ierr = PetscTryMethod(pc,"PCRedundantGetOperators_C",(PC,Mat*,Mat*),(pc,mat,pmat));CHKERRQ(ierr);
4554b9ad928SBarry Smith   PetscFunctionReturn(0);
4564b9ad928SBarry Smith }
4574b9ad928SBarry Smith 
4584b9ad928SBarry Smith /* -------------------------------------------------------------------------------------*/
45937a17b4dSBarry Smith /*MC
46083ab6a24SBarry Smith      PCREDUNDANT - Runs a KSP solver with preconditioner for the entire problem on subgroups of processors
46137a17b4dSBarry Smith 
46283ab6a24SBarry Smith      Options for the redundant preconditioners can be set with -redundant_pc_xxx for the redundant KSP with -redundant_ksp_xxx
46337a17b4dSBarry Smith 
46409391456SBarry Smith   Options Database:
4659b21d695SBarry Smith .  -pc_redundant_number <n> - number of redundant solves, for example if you are using 64 MPI processes and
4669b21d695SBarry Smith                               use an n of 4 there will be 4 parallel solves each on 16 = 64/4 processes.
46709391456SBarry Smith 
46837a17b4dSBarry Smith    Level: intermediate
46937a17b4dSBarry Smith 
47083ab6a24SBarry Smith    Notes: The default KSP is preonly and the default PC is LU.
47183ab6a24SBarry Smith 
47283ab6a24SBarry Smith    Developer Notes: Note that PCSetInitialGuessNonzero()  is not used by this class but likely should be.
4739cfaa89bSBarry Smith 
47437a17b4dSBarry Smith .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PCRedundantSetScatter(),
47583ab6a24SBarry Smith            PCRedundantGetKSP(), PCRedundantGetOperators(), PCRedundantSetNumber()
47637a17b4dSBarry Smith M*/
47737a17b4dSBarry Smith 
4784b9ad928SBarry Smith #undef __FUNCT__
4794b9ad928SBarry Smith #define __FUNCT__ "PCCreate_Redundant"
4808cc058d9SJed Brown PETSC_EXTERN PetscErrorCode PCCreate_Redundant(PC pc)
4814b9ad928SBarry Smith {
482dfbe8321SBarry Smith   PetscErrorCode ierr;
4834b9ad928SBarry Smith   PC_Redundant   *red;
48469db28dcSHong Zhang   PetscMPIInt    size;
4853f457be1SHong Zhang 
4864b9ad928SBarry Smith   PetscFunctionBegin;
48738f2d2fdSLisandro Dalcin   ierr = PetscNewLog(pc,PC_Redundant,&red);CHKERRQ(ierr);
488ce94432eSBarry Smith   ierr = MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);CHKERRQ(ierr);
4892fa5cd67SKarl Rupp 
49069db28dcSHong Zhang   red->nsubcomm       = size;
4914b9ad928SBarry Smith   red->useparallelmat = PETSC_TRUE;
4921fbd8f88SHong Zhang   pc->data            = (void*)red;
4934b9ad928SBarry Smith 
4944b9ad928SBarry Smith   pc->ops->apply          = PCApply_Redundant;
4954b9ad928SBarry Smith   pc->ops->applytranspose = 0;
4964b9ad928SBarry Smith   pc->ops->setup          = PCSetUp_Redundant;
4974b9ad928SBarry Smith   pc->ops->destroy        = PCDestroy_Redundant;
4981ea5a559SBarry Smith   pc->ops->reset          = PCReset_Redundant;
4994b9ad928SBarry Smith   pc->ops->setfromoptions = PCSetFromOptions_Redundant;
5004b9ad928SBarry Smith   pc->ops->view           = PCView_Redundant;
5012fa5cd67SKarl Rupp 
502bdf89e91SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCRedundantSetScatter_C",PCRedundantSetScatter_Redundant);CHKERRQ(ierr);
503bdf89e91SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCRedundantSetNumber_C",PCRedundantSetNumber_Redundant);CHKERRQ(ierr);
504bdf89e91SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCRedundantGetKSP_C",PCRedundantGetKSP_Redundant);CHKERRQ(ierr);
505bdf89e91SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCRedundantGetOperators_C",PCRedundantGetOperators_Redundant);CHKERRQ(ierr);
5064b9ad928SBarry Smith   PetscFunctionReturn(0);
5074b9ad928SBarry Smith }
508b2573a8aSBarry Smith 
509