xref: /petsc/src/ksp/pc/impls/redundant/redundant.c (revision c1619fb62554bd9e7cce5c1feff68159f1df823f)
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 */
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);
373f08860eSBarry Smith       ierr = PetscViewerGetSubViewer(viewer,((PetscObject)red->pc)->comm,&subviewer);CHKERRQ(ierr);
38f5dd71faSBarry Smith       if (!red->psubcomm->color) { /* only view first redundant pc */
391575c14dSBarry Smith         ierr = PetscViewerASCIIPushTab(subviewer);CHKERRQ(ierr);
403e065800SHong Zhang         ierr = KSPView(red->ksp,subviewer);CHKERRQ(ierr);
411575c14dSBarry Smith         ierr = PetscViewerASCIIPopTab(subviewer);CHKERRQ(ierr);
424b9ad928SBarry Smith       }
433f08860eSBarry Smith       ierr = PetscViewerRestoreSubViewer(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 
514b9ad928SBarry Smith #undef __FUNCT__
524b9ad928SBarry Smith #define __FUNCT__ "PCSetUp_Redundant"
536849ba73SBarry Smith static PetscErrorCode PCSetUp_Redundant(PC pc)
544b9ad928SBarry Smith {
554b9ad928SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
56dfbe8321SBarry Smith   PetscErrorCode ierr;
571b81debcSHong Zhang   PetscInt       mstart,mend,mlocal,M;
5813f74950SBarry Smith   PetscMPIInt    size;
59ce94432eSBarry Smith   MPI_Comm       comm,subcomm;
60ddc54837SHong Zhang   Vec            x;
613f457be1SHong Zhang 
624b9ad928SBarry Smith   PetscFunctionBegin;
63ce94432eSBarry Smith   ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
64ddc54837SHong Zhang 
65ddc54837SHong Zhang   /* if pmatrix set by user is sequential then we do not need to gather the parallel matrix */
66ddc54837SHong Zhang   ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
67ddc54837SHong Zhang   if (size == 1) red->useparallelmat = PETSC_FALSE;
681fbd8f88SHong Zhang 
694b9ad928SBarry Smith   if (!pc->setupcalled) {
701b81debcSHong Zhang     PetscInt mloc_sub;
7175024027SHong Zhang     if (!red->psubcomm) { /* create red->psubcomm, new ksp and pc over subcomm */
7275024027SHong Zhang       KSP ksp;
7375024027SHong Zhang       ierr = PCRedundantGetKSP(pc,&ksp);CHKERRQ(ierr);
741b81debcSHong Zhang     }
7575024027SHong Zhang     subcomm = PetscSubcommChild(red->psubcomm);
761fbd8f88SHong Zhang 
771b81debcSHong Zhang     if (red->useparallelmat) {
781b81debcSHong Zhang       /* grab the parallel matrix and put it into processors of a subcomminicator */
7953cd1579SHong Zhang       ierr = MatCreateRedundantMatrix(pc->pmat,red->psubcomm->n,subcomm,MAT_INITIAL_MATRIX,&red->pmats);CHKERRQ(ierr);
80b85f2e9bSHong Zhang 
81b85f2e9bSHong Zhang       ierr = MPI_Comm_size(subcomm,&size);CHKERRQ(ierr);
82b85f2e9bSHong Zhang       if (size > 1) {
83b85f2e9bSHong Zhang         PetscBool foundpack;
84b85f2e9bSHong Zhang         ierr = MatGetFactorAvailable(red->pmats,NULL,MAT_FACTOR_LU,&foundpack);CHKERRQ(ierr);
85b85f2e9bSHong Zhang         if (!foundpack) { /* reset default ksp and pc */
86b85f2e9bSHong Zhang           ierr = KSPSetType(red->ksp,KSPGMRES);CHKERRQ(ierr);
87b85f2e9bSHong Zhang           ierr = PCSetType(red->pc,PCBJACOBI);CHKERRQ(ierr);
88*c1619fb6SBarry Smith         } else {
89*c1619fb6SBarry Smith           ierr = PCFactorSetMatSolverPackage(red->pc,NULL);CHKERRQ(ierr);
90b85f2e9bSHong Zhang         }
91b85f2e9bSHong Zhang       }
92b85f2e9bSHong Zhang 
9323ee1639SBarry Smith       ierr = KSPSetOperators(red->ksp,red->pmats,red->pmats);CHKERRQ(ierr);
944b9ad928SBarry Smith 
951b81debcSHong Zhang       /* get working vectors xsub and ysub */
962a7a6963SBarry Smith       ierr = MatCreateVecs(red->pmats,&red->xsub,&red->ysub);CHKERRQ(ierr);
972fa5cd67SKarl Rupp 
988b96b0d2SHong Zhang       /* create working vectors xdup and ydup.
998b96b0d2SHong Zhang        xdup concatenates all xsub's contigously to form a mpi vector over dupcomm  (see PetscSubcommCreate_interlaced())
1008b96b0d2SHong Zhang        ydup concatenates all ysub and has empty local arrays because ysub's arrays will be place into it.
101ce94432eSBarry Smith        Note: we use communicator dupcomm, not PetscObjectComm((PetscObject)pc)! */
1021b81debcSHong Zhang       ierr = MatGetLocalSize(red->pmats,&mloc_sub,NULL);CHKERRQ(ierr);
10336be1a5eSBarry Smith       ierr = VecCreateMPI(PetscSubcommContiguousParent(red->psubcomm),mloc_sub,PETSC_DECIDE,&red->xdup);CHKERRQ(ierr);
10436be1a5eSBarry Smith       ierr = VecCreateMPIWithArray(PetscSubcommContiguousParent(red->psubcomm),1,mloc_sub,PETSC_DECIDE,NULL,&red->ydup);CHKERRQ(ierr);
1053f457be1SHong Zhang 
106f68be91cSHong Zhang       /* create vecscatters */
107f68be91cSHong Zhang       if (!red->scatterin) { /* efficiency of scatterin is independent from psubcomm_type! */
1083f457be1SHong Zhang         IS       is1,is2;
1093f457be1SHong Zhang         PetscInt *idx1,*idx2,i,j,k;
11045fc02eaSBarry Smith 
1112a7a6963SBarry Smith         ierr = MatCreateVecs(pc->pmat,&x,0);CHKERRQ(ierr);
1121b81debcSHong Zhang         ierr = VecGetSize(x,&M);CHKERRQ(ierr);
1131b81debcSHong Zhang         ierr = VecGetOwnershipRange(x,&mstart,&mend);CHKERRQ(ierr);
1141b81debcSHong Zhang         mlocal = mend - mstart;
115dcca6d9dSJed Brown         ierr = PetscMalloc2(red->psubcomm->n*mlocal,&idx1,red->psubcomm->n*mlocal,&idx2);CHKERRQ(ierr);
1163f457be1SHong Zhang         j    = 0;
1171fbd8f88SHong Zhang         for (k=0; k<red->psubcomm->n; k++) {
1183f457be1SHong Zhang           for (i=mstart; i<mend; i++) {
1193f457be1SHong Zhang             idx1[j]   = i;
120ddc54837SHong Zhang             idx2[j++] = i + M*k;
1213f457be1SHong Zhang           }
1223f457be1SHong Zhang         }
12370b3c8c7SBarry Smith         ierr = ISCreateGeneral(comm,red->psubcomm->n*mlocal,idx1,PETSC_COPY_VALUES,&is1);CHKERRQ(ierr);
12470b3c8c7SBarry Smith         ierr = ISCreateGeneral(comm,red->psubcomm->n*mlocal,idx2,PETSC_COPY_VALUES,&is2);CHKERRQ(ierr);
125ddc54837SHong Zhang         ierr = VecScatterCreate(x,is1,red->xdup,is2,&red->scatterin);CHKERRQ(ierr);
126fcfd50ebSBarry Smith         ierr = ISDestroy(&is1);CHKERRQ(ierr);
127fcfd50ebSBarry Smith         ierr = ISDestroy(&is2);CHKERRQ(ierr);
1283f457be1SHong Zhang 
1296909748bSHong Zhang         /* Impl below is good for PETSC_SUBCOMM_INTERLACED (no inter-process communication) and PETSC_SUBCOMM_CONTIGUOUS (communication within subcomm) */
130ddc54837SHong Zhang         ierr = ISCreateStride(comm,mlocal,mstart+ red->psubcomm->color*M,1,&is1);CHKERRQ(ierr);
1313f457be1SHong Zhang         ierr = ISCreateStride(comm,mlocal,mstart,1,&is2);CHKERRQ(ierr);
132ddc54837SHong Zhang         ierr = VecScatterCreate(red->xdup,is1,x,is2,&red->scatterout);CHKERRQ(ierr);
133fcfd50ebSBarry Smith         ierr = ISDestroy(&is1);CHKERRQ(ierr);
134fcfd50ebSBarry Smith         ierr = ISDestroy(&is2);CHKERRQ(ierr);
1351d79065fSBarry Smith         ierr = PetscFree2(idx1,idx2);CHKERRQ(ierr);
136ddc54837SHong Zhang         ierr = VecDestroy(&x);CHKERRQ(ierr);
1371b81debcSHong Zhang       }
138ab661555SHong Zhang     } else { /* !red->useparallelmat */
13923ee1639SBarry Smith       ierr = KSPSetOperators(red->ksp,pc->mat,pc->pmat);CHKERRQ(ierr);
1401b81debcSHong Zhang     }
141ab661555SHong Zhang   } else { /* pc->setupcalled */
1424b9ad928SBarry Smith     if (red->useparallelmat) {
143ab661555SHong Zhang       MatReuse       reuse;
1441b81debcSHong Zhang       /* grab the parallel matrix and put it into processors of a subcomminicator */
1451b81debcSHong Zhang       /*--------------------------------------------------------------------------*/
146ab661555SHong Zhang       if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1474b9ad928SBarry Smith         /* destroy old matrices */
1486bf464f9SBarry Smith         ierr  = MatDestroy(&red->pmats);CHKERRQ(ierr);
149ab661555SHong Zhang         reuse = MAT_INITIAL_MATRIX;
1504b9ad928SBarry Smith       } else {
151ab661555SHong Zhang         reuse = MAT_REUSE_MATRIX;
152ab661555SHong Zhang       }
153306c2d5bSBarry Smith       ierr = MatCreateRedundantMatrix(pc->pmat,red->psubcomm->n,PetscSubcommChild(red->psubcomm),reuse,&red->pmats);CHKERRQ(ierr);
15423ee1639SBarry Smith       ierr = KSPSetOperators(red->ksp,red->pmats,red->pmats);CHKERRQ(ierr);
155ab661555SHong Zhang     } else { /* !red->useparallelmat */
15623ee1639SBarry Smith       ierr = KSPSetOperators(red->ksp,pc->mat,pc->pmat);CHKERRQ(ierr);
1574b9ad928SBarry Smith     }
158ab661555SHong Zhang   }
1591b81debcSHong Zhang 
1600c24e6a1SHong Zhang   if (pc->setfromoptionscalled) {
1613e065800SHong Zhang     ierr = KSPSetFromOptions(red->ksp);CHKERRQ(ierr);
1620c24e6a1SHong Zhang   }
1633e065800SHong Zhang   ierr = KSPSetUp(red->ksp);CHKERRQ(ierr);
1644b9ad928SBarry Smith   PetscFunctionReturn(0);
1654b9ad928SBarry Smith }
1664b9ad928SBarry Smith 
1674b9ad928SBarry Smith #undef __FUNCT__
1684b9ad928SBarry Smith #define __FUNCT__ "PCApply_Redundant"
1696849ba73SBarry Smith static PetscErrorCode PCApply_Redundant(PC pc,Vec x,Vec y)
1704b9ad928SBarry Smith {
1714b9ad928SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
172dfbe8321SBarry Smith   PetscErrorCode ierr;
1733f457be1SHong Zhang   PetscScalar    *array;
1744b9ad928SBarry Smith 
1754b9ad928SBarry Smith   PetscFunctionBegin;
176ddc54837SHong Zhang   if (!red->useparallelmat) {
177ddc54837SHong Zhang     ierr = KSPSolve(red->ksp,x,y);CHKERRQ(ierr);
178ddc54837SHong Zhang     PetscFunctionReturn(0);
179ddc54837SHong Zhang   }
180ddc54837SHong Zhang 
1813f457be1SHong Zhang   /* scatter x to xdup */
182ca9f406cSSatish Balay   ierr = VecScatterBegin(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
183ca9f406cSSatish Balay   ierr = VecScatterEnd(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1843f457be1SHong Zhang 
1853f457be1SHong Zhang   /* place xdup's local array into xsub */
1863f457be1SHong Zhang   ierr = VecGetArray(red->xdup,&array);CHKERRQ(ierr);
1873f457be1SHong Zhang   ierr = VecPlaceArray(red->xsub,(const PetscScalar*)array);CHKERRQ(ierr);
1884b9ad928SBarry Smith 
1894b9ad928SBarry Smith   /* apply preconditioner on each processor */
19083ab6a24SBarry Smith   ierr = KSPSolve(red->ksp,red->xsub,red->ysub);CHKERRQ(ierr);
1913f457be1SHong Zhang   ierr = VecResetArray(red->xsub);CHKERRQ(ierr);
1923f457be1SHong Zhang   ierr = VecRestoreArray(red->xdup,&array);CHKERRQ(ierr);
1934b9ad928SBarry Smith 
1943f457be1SHong Zhang   /* place ysub's local array into ydup */
1953f457be1SHong Zhang   ierr = VecGetArray(red->ysub,&array);CHKERRQ(ierr);
1963f457be1SHong Zhang   ierr = VecPlaceArray(red->ydup,(const PetscScalar*)array);CHKERRQ(ierr);
1973f457be1SHong Zhang 
1983f457be1SHong Zhang   /* scatter ydup to y */
199ca9f406cSSatish Balay   ierr = VecScatterBegin(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
200ca9f406cSSatish Balay   ierr = VecScatterEnd(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
2013f457be1SHong Zhang   ierr = VecResetArray(red->ydup);CHKERRQ(ierr);
2023f457be1SHong Zhang   ierr = VecRestoreArray(red->ysub,&array);CHKERRQ(ierr);
2034b9ad928SBarry Smith   PetscFunctionReturn(0);
2044b9ad928SBarry Smith }
2054b9ad928SBarry Smith 
2064b9ad928SBarry Smith #undef __FUNCT__
207d88bfacbSStefano Zampini #define __FUNCT__ "PCApplyTranspose_Redundant"
208d88bfacbSStefano Zampini static PetscErrorCode PCApplyTranspose_Redundant(PC pc,Vec x,Vec y)
209d88bfacbSStefano Zampini {
210d88bfacbSStefano Zampini   PC_Redundant   *red = (PC_Redundant*)pc->data;
211d88bfacbSStefano Zampini   PetscErrorCode ierr;
212d88bfacbSStefano Zampini   PetscScalar    *array;
213d88bfacbSStefano Zampini 
214d88bfacbSStefano Zampini   PetscFunctionBegin;
215d88bfacbSStefano Zampini   if (!red->useparallelmat) {
216d88bfacbSStefano Zampini     ierr = KSPSolveTranspose(red->ksp,x,y);CHKERRQ(ierr);
217d88bfacbSStefano Zampini     PetscFunctionReturn(0);
218d88bfacbSStefano Zampini   }
219d88bfacbSStefano Zampini 
220d88bfacbSStefano Zampini   /* scatter x to xdup */
221d88bfacbSStefano Zampini   ierr = VecScatterBegin(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
222d88bfacbSStefano Zampini   ierr = VecScatterEnd(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
223d88bfacbSStefano Zampini 
224d88bfacbSStefano Zampini   /* place xdup's local array into xsub */
225d88bfacbSStefano Zampini   ierr = VecGetArray(red->xdup,&array);CHKERRQ(ierr);
226d88bfacbSStefano Zampini   ierr = VecPlaceArray(red->xsub,(const PetscScalar*)array);CHKERRQ(ierr);
227d88bfacbSStefano Zampini 
228d88bfacbSStefano Zampini   /* apply preconditioner on each processor */
229d88bfacbSStefano Zampini   ierr = KSPSolveTranspose(red->ksp,red->xsub,red->ysub);CHKERRQ(ierr);
230d88bfacbSStefano Zampini   ierr = VecResetArray(red->xsub);CHKERRQ(ierr);
231d88bfacbSStefano Zampini   ierr = VecRestoreArray(red->xdup,&array);CHKERRQ(ierr);
232d88bfacbSStefano Zampini 
233d88bfacbSStefano Zampini   /* place ysub's local array into ydup */
234d88bfacbSStefano Zampini   ierr = VecGetArray(red->ysub,&array);CHKERRQ(ierr);
235d88bfacbSStefano Zampini   ierr = VecPlaceArray(red->ydup,(const PetscScalar*)array);CHKERRQ(ierr);
236d88bfacbSStefano Zampini 
237d88bfacbSStefano Zampini   /* scatter ydup to y */
238d88bfacbSStefano Zampini   ierr = VecScatterBegin(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
239d88bfacbSStefano Zampini   ierr = VecScatterEnd(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
240d88bfacbSStefano Zampini   ierr = VecResetArray(red->ydup);CHKERRQ(ierr);
241d88bfacbSStefano Zampini   ierr = VecRestoreArray(red->ysub,&array);CHKERRQ(ierr);
242d88bfacbSStefano Zampini   PetscFunctionReturn(0);
243d88bfacbSStefano Zampini }
244d88bfacbSStefano Zampini 
245d88bfacbSStefano Zampini #undef __FUNCT__
2461ea5a559SBarry Smith #define __FUNCT__ "PCReset_Redundant"
2471ea5a559SBarry Smith static PetscErrorCode PCReset_Redundant(PC pc)
2484b9ad928SBarry Smith {
2494b9ad928SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
250dfbe8321SBarry Smith   PetscErrorCode ierr;
2514b9ad928SBarry Smith 
2524b9ad928SBarry Smith   PetscFunctionBegin;
2531b81debcSHong Zhang   if (red->useparallelmat) {
2546bf464f9SBarry Smith     ierr = VecScatterDestroy(&red->scatterin);CHKERRQ(ierr);
2556bf464f9SBarry Smith     ierr = VecScatterDestroy(&red->scatterout);CHKERRQ(ierr);
2566bf464f9SBarry Smith     ierr = VecDestroy(&red->ysub);CHKERRQ(ierr);
2576bf464f9SBarry Smith     ierr = VecDestroy(&red->xsub);CHKERRQ(ierr);
2586bf464f9SBarry Smith     ierr = VecDestroy(&red->xdup);CHKERRQ(ierr);
2596bf464f9SBarry Smith     ierr = VecDestroy(&red->ydup);CHKERRQ(ierr);
2601b81debcSHong Zhang   }
2616bf464f9SBarry Smith   ierr = MatDestroy(&red->pmats);CHKERRQ(ierr);
2621b81debcSHong Zhang   ierr = KSPReset(red->ksp);CHKERRQ(ierr);
2631ea5a559SBarry Smith   PetscFunctionReturn(0);
2641ea5a559SBarry Smith }
2651ea5a559SBarry Smith 
2661ea5a559SBarry Smith #undef __FUNCT__
2671ea5a559SBarry Smith #define __FUNCT__ "PCDestroy_Redundant"
2681ea5a559SBarry Smith static PetscErrorCode PCDestroy_Redundant(PC pc)
2691ea5a559SBarry Smith {
2701ea5a559SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
2711ea5a559SBarry Smith   PetscErrorCode ierr;
2721ea5a559SBarry Smith 
2731ea5a559SBarry Smith   PetscFunctionBegin;
2741ea5a559SBarry Smith   ierr = PCReset_Redundant(pc);CHKERRQ(ierr);
2756bf464f9SBarry Smith   ierr = KSPDestroy(&red->ksp);CHKERRQ(ierr);
2766bf464f9SBarry Smith   ierr = PetscSubcommDestroy(&red->psubcomm);CHKERRQ(ierr);
277c31cb41cSBarry Smith   ierr = PetscFree(pc->data);CHKERRQ(ierr);
2784b9ad928SBarry Smith   PetscFunctionReturn(0);
2794b9ad928SBarry Smith }
2804b9ad928SBarry Smith 
2814b9ad928SBarry Smith #undef __FUNCT__
2824b9ad928SBarry Smith #define __FUNCT__ "PCSetFromOptions_Redundant"
2834416b707SBarry Smith static PetscErrorCode PCSetFromOptions_Redundant(PetscOptionItems *PetscOptionsObject,PC pc)
2844b9ad928SBarry Smith {
285a98ce0f4SHong Zhang   PetscErrorCode ierr;
286a98ce0f4SHong Zhang   PC_Redundant   *red = (PC_Redundant*)pc->data;
287a98ce0f4SHong Zhang 
2884b9ad928SBarry Smith   PetscFunctionBegin;
289e55864a3SBarry Smith   ierr = PetscOptionsHead(PetscOptionsObject,"Redundant options");CHKERRQ(ierr);
29009a6bc64SHong Zhang   ierr = PetscOptionsInt("-pc_redundant_number","Number of redundant pc","PCRedundantSetNumber",red->nsubcomm,&red->nsubcomm,0);CHKERRQ(ierr);
291a98ce0f4SHong Zhang   ierr = PetscOptionsTail();CHKERRQ(ierr);
2924b9ad928SBarry Smith   PetscFunctionReturn(0);
2934b9ad928SBarry Smith }
2944b9ad928SBarry Smith 
2954b9ad928SBarry Smith #undef __FUNCT__
29609a6bc64SHong Zhang #define __FUNCT__ "PCRedundantSetNumber_Redundant"
297f7a08781SBarry Smith static PetscErrorCode PCRedundantSetNumber_Redundant(PC pc,PetscInt nreds)
29809a6bc64SHong Zhang {
29909a6bc64SHong Zhang   PC_Redundant *red = (PC_Redundant*)pc->data;
30009a6bc64SHong Zhang 
30109a6bc64SHong Zhang   PetscFunctionBegin;
30209a6bc64SHong Zhang   red->nsubcomm = nreds;
30309a6bc64SHong Zhang   PetscFunctionReturn(0);
30409a6bc64SHong Zhang }
30509a6bc64SHong Zhang 
30609a6bc64SHong Zhang #undef __FUNCT__
30709a6bc64SHong Zhang #define __FUNCT__ "PCRedundantSetNumber"
30809a6bc64SHong Zhang /*@
30909a6bc64SHong Zhang    PCRedundantSetNumber - Sets the number of redundant preconditioner contexts.
31009a6bc64SHong Zhang 
3113f9fe445SBarry Smith    Logically Collective on PC
31209a6bc64SHong Zhang 
31309a6bc64SHong Zhang    Input Parameters:
31409a6bc64SHong Zhang +  pc - the preconditioner context
3159b21d695SBarry Smith -  nredundant - number of redundant preconditioner contexts; for example if you are using 64 MPI processes and
3169b21d695SBarry Smith                               use an nredundant of 4 there will be 4 parallel solves each on 16 = 64/4 processes.
31709a6bc64SHong Zhang 
31809a6bc64SHong Zhang    Level: advanced
31909a6bc64SHong Zhang 
32009a6bc64SHong Zhang .keywords: PC, redundant solve
32109a6bc64SHong Zhang @*/
3227087cfbeSBarry Smith PetscErrorCode PCRedundantSetNumber(PC pc,PetscInt nredundant)
32309a6bc64SHong Zhang {
3244ac538c5SBarry Smith   PetscErrorCode ierr;
32509a6bc64SHong Zhang 
32609a6bc64SHong Zhang   PetscFunctionBegin;
3270700a824SBarry Smith   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
328ce94432eSBarry Smith   if (nredundant <= 0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "num of redundant pc %D must be positive",nredundant);
3294ac538c5SBarry Smith   ierr = PetscTryMethod(pc,"PCRedundantSetNumber_C",(PC,PetscInt),(pc,nredundant));CHKERRQ(ierr);
33009a6bc64SHong Zhang   PetscFunctionReturn(0);
33109a6bc64SHong Zhang }
33209a6bc64SHong Zhang 
33309a6bc64SHong Zhang #undef __FUNCT__
3344b9ad928SBarry Smith #define __FUNCT__ "PCRedundantSetScatter_Redundant"
335f7a08781SBarry Smith static PetscErrorCode PCRedundantSetScatter_Redundant(PC pc,VecScatter in,VecScatter out)
3364b9ad928SBarry Smith {
3374b9ad928SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
338dfbe8321SBarry Smith   PetscErrorCode ierr;
3394b9ad928SBarry Smith 
3404b9ad928SBarry Smith   PetscFunctionBegin;
3414b9ad928SBarry Smith   ierr = PetscObjectReference((PetscObject)in);CHKERRQ(ierr);
3426bf464f9SBarry Smith   ierr = VecScatterDestroy(&red->scatterin);CHKERRQ(ierr);
3432fa5cd67SKarl Rupp 
344c3122656SLisandro Dalcin   red->scatterin  = in;
3452fa5cd67SKarl Rupp 
3464b9ad928SBarry Smith   ierr = PetscObjectReference((PetscObject)out);CHKERRQ(ierr);
3476bf464f9SBarry Smith   ierr = VecScatterDestroy(&red->scatterout);CHKERRQ(ierr);
348c3122656SLisandro Dalcin   red->scatterout = out;
3494b9ad928SBarry Smith   PetscFunctionReturn(0);
3504b9ad928SBarry Smith }
3514b9ad928SBarry Smith 
3524b9ad928SBarry Smith #undef __FUNCT__
3534b9ad928SBarry Smith #define __FUNCT__ "PCRedundantSetScatter"
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 .keywords: PC, redundant solve
3694b9ad928SBarry Smith @*/
3707087cfbeSBarry Smith PetscErrorCode PCRedundantSetScatter(PC pc,VecScatter in,VecScatter out)
3714b9ad928SBarry Smith {
3724ac538c5SBarry Smith   PetscErrorCode ierr;
3734b9ad928SBarry Smith 
3744b9ad928SBarry Smith   PetscFunctionBegin;
3750700a824SBarry Smith   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
3760700a824SBarry Smith   PetscValidHeaderSpecific(in,VEC_SCATTER_CLASSID,2);
3770700a824SBarry Smith   PetscValidHeaderSpecific(out,VEC_SCATTER_CLASSID,3);
3784ac538c5SBarry Smith   ierr = PetscTryMethod(pc,"PCRedundantSetScatter_C",(PC,VecScatter,VecScatter),(pc,in,out));CHKERRQ(ierr);
3794b9ad928SBarry Smith   PetscFunctionReturn(0);
3804b9ad928SBarry Smith }
3814b9ad928SBarry Smith 
3824b9ad928SBarry Smith #undef __FUNCT__
38383ab6a24SBarry Smith #define __FUNCT__ "PCRedundantGetKSP_Redundant"
384f7a08781SBarry Smith static PetscErrorCode PCRedundantGetKSP_Redundant(PC pc,KSP *innerksp)
3854b9ad928SBarry Smith {
3865f06b7aaSBarry Smith   PetscErrorCode ierr;
3874b9ad928SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
38875024027SHong Zhang   MPI_Comm       comm,subcomm;
38975024027SHong Zhang   const char     *prefix;
3904b9ad928SBarry Smith 
3914b9ad928SBarry Smith   PetscFunctionBegin;
39275024027SHong Zhang   if (!red->psubcomm) {
393e5acf8a4SHong Zhang     ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr);
394e5acf8a4SHong Zhang 
39575024027SHong Zhang     ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
39675024027SHong Zhang     ierr = PetscSubcommCreate(comm,&red->psubcomm);CHKERRQ(ierr);
39775024027SHong Zhang     ierr = PetscSubcommSetNumber(red->psubcomm,red->nsubcomm);CHKERRQ(ierr);
39875024027SHong Zhang     ierr = PetscSubcommSetType(red->psubcomm,PETSC_SUBCOMM_CONTIGUOUS);CHKERRQ(ierr);
399e5acf8a4SHong Zhang 
400e5acf8a4SHong Zhang     ierr = PetscSubcommSetOptionsPrefix(red->psubcomm,prefix);CHKERRQ(ierr);
401e5acf8a4SHong Zhang     ierr = PetscSubcommSetFromOptions(red->psubcomm);CHKERRQ(ierr);
40275024027SHong Zhang     ierr = PetscLogObjectMemory((PetscObject)pc,sizeof(PetscSubcomm));CHKERRQ(ierr);
40375024027SHong Zhang 
40475024027SHong Zhang     /* create a new PC that processors in each subcomm have copy of */
40575024027SHong Zhang     subcomm = PetscSubcommChild(red->psubcomm);
40675024027SHong Zhang 
40775024027SHong Zhang     ierr = KSPCreate(subcomm,&red->ksp);CHKERRQ(ierr);
40875024027SHong Zhang     ierr = KSPSetErrorIfNotConverged(red->ksp,pc->erroriffailure);CHKERRQ(ierr);
40975024027SHong Zhang     ierr = PetscObjectIncrementTabLevel((PetscObject)red->ksp,(PetscObject)pc,1);CHKERRQ(ierr);
41075024027SHong Zhang     ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)red->ksp);CHKERRQ(ierr);
41175024027SHong Zhang     ierr = KSPSetType(red->ksp,KSPPREONLY);CHKERRQ(ierr);
41275024027SHong Zhang     ierr = KSPGetPC(red->ksp,&red->pc);CHKERRQ(ierr);
41375024027SHong Zhang     ierr = PCSetType(red->pc,PCLU);CHKERRQ(ierr);
41475024027SHong Zhang 
41575024027SHong Zhang     ierr = KSPSetOptionsPrefix(red->ksp,prefix);CHKERRQ(ierr);
41675024027SHong Zhang     ierr = KSPAppendOptionsPrefix(red->ksp,"redundant_");CHKERRQ(ierr);
41775024027SHong Zhang   }
41883ab6a24SBarry Smith   *innerksp = red->ksp;
4194b9ad928SBarry Smith   PetscFunctionReturn(0);
4204b9ad928SBarry Smith }
4214b9ad928SBarry Smith 
4224b9ad928SBarry Smith #undef __FUNCT__
42383ab6a24SBarry Smith #define __FUNCT__ "PCRedundantGetKSP"
4244b9ad928SBarry Smith /*@
42583ab6a24SBarry Smith    PCRedundantGetKSP - Gets the less parallel KSP created by the redundant PC.
4264b9ad928SBarry Smith 
4274b9ad928SBarry Smith    Not Collective
4284b9ad928SBarry Smith 
4294b9ad928SBarry Smith    Input Parameter:
4304b9ad928SBarry Smith .  pc - the preconditioner context
4314b9ad928SBarry Smith 
4324b9ad928SBarry Smith    Output Parameter:
43383ab6a24SBarry Smith .  innerksp - the KSP on the smaller set of processes
4344b9ad928SBarry Smith 
4354b9ad928SBarry Smith    Level: advanced
4364b9ad928SBarry Smith 
4374b9ad928SBarry Smith .keywords: PC, redundant solve
4384b9ad928SBarry Smith @*/
43983ab6a24SBarry Smith PetscErrorCode PCRedundantGetKSP(PC pc,KSP *innerksp)
4404b9ad928SBarry Smith {
4414ac538c5SBarry Smith   PetscErrorCode ierr;
4424b9ad928SBarry Smith 
4434b9ad928SBarry Smith   PetscFunctionBegin;
4440700a824SBarry Smith   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
44583ab6a24SBarry Smith   PetscValidPointer(innerksp,2);
44683ab6a24SBarry Smith   ierr = PetscTryMethod(pc,"PCRedundantGetKSP_C",(PC,KSP*),(pc,innerksp));CHKERRQ(ierr);
4474b9ad928SBarry Smith   PetscFunctionReturn(0);
4484b9ad928SBarry Smith }
4494b9ad928SBarry Smith 
4504b9ad928SBarry Smith #undef __FUNCT__
4514b9ad928SBarry Smith #define __FUNCT__ "PCRedundantGetOperators_Redundant"
452f7a08781SBarry Smith static PetscErrorCode PCRedundantGetOperators_Redundant(PC pc,Mat *mat,Mat *pmat)
4534b9ad928SBarry Smith {
4544b9ad928SBarry Smith   PC_Redundant *red = (PC_Redundant*)pc->data;
4554b9ad928SBarry Smith 
4564b9ad928SBarry Smith   PetscFunctionBegin;
457b3804887SHong Zhang   if (mat)  *mat  = red->pmats;
458b3804887SHong Zhang   if (pmat) *pmat = red->pmats;
4594b9ad928SBarry Smith   PetscFunctionReturn(0);
4604b9ad928SBarry Smith }
4614b9ad928SBarry Smith 
4624b9ad928SBarry Smith #undef __FUNCT__
4634b9ad928SBarry Smith #define __FUNCT__ "PCRedundantGetOperators"
4644b9ad928SBarry Smith /*@
4654b9ad928SBarry Smith    PCRedundantGetOperators - gets the sequential matrix and preconditioner matrix
4664b9ad928SBarry Smith 
4674b9ad928SBarry Smith    Not Collective
4684b9ad928SBarry Smith 
4694b9ad928SBarry Smith    Input Parameter:
4704b9ad928SBarry Smith .  pc - the preconditioner context
4714b9ad928SBarry Smith 
4724b9ad928SBarry Smith    Output Parameters:
4734b9ad928SBarry Smith +  mat - the matrix
4744b9ad928SBarry Smith -  pmat - the (possibly different) preconditioner matrix
4754b9ad928SBarry Smith 
4764b9ad928SBarry Smith    Level: advanced
4774b9ad928SBarry Smith 
4784b9ad928SBarry Smith .keywords: PC, redundant solve
4794b9ad928SBarry Smith @*/
4807087cfbeSBarry Smith PetscErrorCode PCRedundantGetOperators(PC pc,Mat *mat,Mat *pmat)
4814b9ad928SBarry Smith {
4824ac538c5SBarry Smith   PetscErrorCode ierr;
4834b9ad928SBarry Smith 
4844b9ad928SBarry Smith   PetscFunctionBegin;
4850700a824SBarry Smith   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
4864482741eSBarry Smith   if (mat)  PetscValidPointer(mat,2);
4874482741eSBarry Smith   if (pmat) PetscValidPointer(pmat,3);
4884ac538c5SBarry Smith   ierr = PetscTryMethod(pc,"PCRedundantGetOperators_C",(PC,Mat*,Mat*),(pc,mat,pmat));CHKERRQ(ierr);
4894b9ad928SBarry Smith   PetscFunctionReturn(0);
4904b9ad928SBarry Smith }
4914b9ad928SBarry Smith 
4924b9ad928SBarry Smith /* -------------------------------------------------------------------------------------*/
49337a17b4dSBarry Smith /*MC
49483ab6a24SBarry Smith      PCREDUNDANT - Runs a KSP solver with preconditioner for the entire problem on subgroups of processors
49537a17b4dSBarry Smith 
49683ab6a24SBarry Smith      Options for the redundant preconditioners can be set with -redundant_pc_xxx for the redundant KSP with -redundant_ksp_xxx
49737a17b4dSBarry Smith 
49809391456SBarry Smith   Options Database:
4999b21d695SBarry Smith .  -pc_redundant_number <n> - number of redundant solves, for example if you are using 64 MPI processes and
5009b21d695SBarry Smith                               use an n of 4 there will be 4 parallel solves each on 16 = 64/4 processes.
50109391456SBarry Smith 
50237a17b4dSBarry Smith    Level: intermediate
50337a17b4dSBarry Smith 
50483ab6a24SBarry Smith    Notes: The default KSP is preonly and the default PC is LU.
50583ab6a24SBarry Smith 
50683ab6a24SBarry Smith    Developer Notes: Note that PCSetInitialGuessNonzero()  is not used by this class but likely should be.
5079cfaa89bSBarry Smith 
50837a17b4dSBarry Smith .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PCRedundantSetScatter(),
50983ab6a24SBarry Smith            PCRedundantGetKSP(), PCRedundantGetOperators(), PCRedundantSetNumber()
51037a17b4dSBarry Smith M*/
51137a17b4dSBarry Smith 
5124b9ad928SBarry Smith #undef __FUNCT__
5134b9ad928SBarry Smith #define __FUNCT__ "PCCreate_Redundant"
5148cc058d9SJed Brown PETSC_EXTERN PetscErrorCode PCCreate_Redundant(PC pc)
5154b9ad928SBarry Smith {
516dfbe8321SBarry Smith   PetscErrorCode ierr;
5174b9ad928SBarry Smith   PC_Redundant   *red;
51869db28dcSHong Zhang   PetscMPIInt    size;
5193f457be1SHong Zhang 
5204b9ad928SBarry Smith   PetscFunctionBegin;
521b00a9115SJed Brown   ierr = PetscNewLog(pc,&red);CHKERRQ(ierr);
522ce94432eSBarry Smith   ierr = MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);CHKERRQ(ierr);
5232fa5cd67SKarl Rupp 
52469db28dcSHong Zhang   red->nsubcomm       = size;
5254b9ad928SBarry Smith   red->useparallelmat = PETSC_TRUE;
5261fbd8f88SHong Zhang   pc->data            = (void*)red;
5274b9ad928SBarry Smith 
5284b9ad928SBarry Smith   pc->ops->apply          = PCApply_Redundant;
529d88bfacbSStefano Zampini   pc->ops->applytranspose = PCApplyTranspose_Redundant;
5304b9ad928SBarry Smith   pc->ops->setup          = PCSetUp_Redundant;
5314b9ad928SBarry Smith   pc->ops->destroy        = PCDestroy_Redundant;
5321ea5a559SBarry Smith   pc->ops->reset          = PCReset_Redundant;
5334b9ad928SBarry Smith   pc->ops->setfromoptions = PCSetFromOptions_Redundant;
5344b9ad928SBarry Smith   pc->ops->view           = PCView_Redundant;
5352fa5cd67SKarl Rupp 
536bdf89e91SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCRedundantSetScatter_C",PCRedundantSetScatter_Redundant);CHKERRQ(ierr);
537bdf89e91SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCRedundantSetNumber_C",PCRedundantSetNumber_Redundant);CHKERRQ(ierr);
538bdf89e91SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCRedundantGetKSP_C",PCRedundantGetKSP_Redundant);CHKERRQ(ierr);
539bdf89e91SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCRedundantGetOperators_C",PCRedundantGetOperators_Redundant);CHKERRQ(ierr);
5404b9ad928SBarry Smith   PetscFunctionReturn(0);
5414b9ad928SBarry Smith }
542b2573a8aSBarry Smith 
543