xref: /petsc/src/ksp/pc/impls/redundant/redundant.c (revision 69db28dc0cb02ea7f128a31d92966470b5957e3f)
1dba47a55SKris Buschelman #define PETSCKSP_DLL
2dba47a55SKris Buschelman 
34b9ad928SBarry Smith /*
43f457be1SHong Zhang   This file defines a "solve the problem redundantly on each subgroup of processor" preconditioner.
54b9ad928SBarry Smith */
66356e834SBarry Smith #include "private/pcimpl.h"     /*I "petscpc.h" I*/
74b9ad928SBarry Smith #include "petscksp.h"
84b9ad928SBarry Smith 
94b9ad928SBarry Smith typedef struct {
101fbd8f88SHong Zhang   MPI_Comm   parent;      /* parent communicator */
111fbd8f88SHong Zhang   MPI_Comm   dupparent;   /* duplicate parent communicator, under which the processors of this subcomm have contiguous rank */
121fbd8f88SHong Zhang   MPI_Comm   comm;        /* this communicator */
131fbd8f88SHong Zhang   PetscInt   n;           /* num of subcommunicators under the parent communicator */
141fbd8f88SHong Zhang   PetscInt   color;       /* color of processors belong to this communicator */
151fbd8f88SHong Zhang } PetscSubcomm;
161fbd8f88SHong Zhang 
171fbd8f88SHong Zhang #undef __FUNCT__
181fbd8f88SHong Zhang #define __FUNCT__ "PetscSubcommDestroy"
191fbd8f88SHong Zhang PetscErrorCode PetscSubcommDestroy(PetscSubcomm *psubcomm)
201fbd8f88SHong Zhang {
211fbd8f88SHong Zhang   PetscErrorCode ierr;
221fbd8f88SHong Zhang 
231fbd8f88SHong Zhang   PetscFunctionBegin;
241fbd8f88SHong Zhang   ierr = PetscFree(psubcomm);CHKERRQ(ierr);
251fbd8f88SHong Zhang   PetscFunctionReturn(0);
261fbd8f88SHong Zhang }
271fbd8f88SHong Zhang 
281fbd8f88SHong Zhang /*--------------------------------------------------------------------------------------------------
291fbd8f88SHong Zhang     To avoid data scattering from subcomm back to original comm, we create subcomm by iteratively taking a
301fbd8f88SHong Zhang      processe into a subcomm.
311fbd8f88SHong Zhang      An example: size=4, nsubcomm=3
321fbd8f88SHong Zhang       pc->comm:
331fbd8f88SHong Zhang       rank:     [0]  [1]  [2]  [3]
341fbd8f88SHong Zhang       color:     0    1    2    0
351fbd8f88SHong Zhang 
361fbd8f88SHong Zhang      subcomm:
371fbd8f88SHong Zhang       subrank:  [0]  [0]  [0]  [1]
381fbd8f88SHong Zhang 
391fbd8f88SHong Zhang      dupcomm:
401fbd8f88SHong Zhang       duprank:  [0]  [2]  [3]  [1]
411fbd8f88SHong Zhang 
421fbd8f88SHong Zhang      Here, subcomm[color = 0] has subsize=2, owns process [0] and [3]
431fbd8f88SHong Zhang            subcomm[color = 1] has subsize=1, owns process [1]
441fbd8f88SHong Zhang            subcomm[color = 2] has subsize=1, owns process [2]
451fbd8f88SHong Zhang           dupcomm has same number of processes as pc->comm, and its duprank maps
461fbd8f88SHong Zhang           processes in subcomm contiguously into a 1d array:
471fbd8f88SHong Zhang             duprank: [0] [1]      [2]         [3]
481fbd8f88SHong Zhang             rank:    [0] [3]      [1]         [2]
491fbd8f88SHong Zhang                     subcomm[0] subcomm[1]  subcomm[2]
501fbd8f88SHong Zhang      ----------------------------------------------------------------------------------------*/
511fbd8f88SHong Zhang #undef __FUNCT__
521fbd8f88SHong Zhang #define __FUNCT__ "PetscSubcommCreate"
531fbd8f88SHong Zhang PetscErrorCode PetscSubcommCreate(MPI_Comm comm,PetscInt nsubcomm,PetscSubcomm **psubcomm)
541fbd8f88SHong Zhang {
551fbd8f88SHong Zhang   PetscErrorCode ierr;
561fbd8f88SHong Zhang   PetscMPIInt    rank,size,*subsize,duprank,subrank;
571fbd8f88SHong Zhang   PetscInt       np_subcomm,nleftover,i,j,color;
581fbd8f88SHong Zhang   MPI_Comm       subcomm=0,dupcomm=0;
591fbd8f88SHong Zhang   const char     *prefix;
601fbd8f88SHong Zhang   PetscSubcomm   *psubcomm_tmp;
611fbd8f88SHong Zhang 
621fbd8f88SHong Zhang   PetscFunctionBegin;
631fbd8f88SHong Zhang   ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
641fbd8f88SHong Zhang   ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
651fbd8f88SHong Zhang   if (nsubcomm < 1 || nsubcomm > size) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Num of subcommunicators %D cannot be < 1 or > input comm size %D",nsubcomm,size);
661fbd8f88SHong Zhang 
671fbd8f88SHong Zhang   /* get size of each subcommunicator */
681fbd8f88SHong Zhang   ierr = PetscMalloc((1+nsubcomm)*sizeof(PetscMPIInt),&subsize);CHKERRQ(ierr);
691fbd8f88SHong Zhang   np_subcomm = size/nsubcomm;
701fbd8f88SHong Zhang   nleftover  = size - nsubcomm*np_subcomm;
711fbd8f88SHong Zhang   for (i=0; i<nsubcomm; i++){
721fbd8f88SHong Zhang     subsize[i] = np_subcomm;
731fbd8f88SHong Zhang     if (i<nleftover) subsize[i]++;
741fbd8f88SHong Zhang   }
751fbd8f88SHong Zhang 
761fbd8f88SHong Zhang   /* find color for this proc */
771fbd8f88SHong Zhang   color   = rank%nsubcomm;
781fbd8f88SHong Zhang   subrank = rank/nsubcomm;
791fbd8f88SHong Zhang 
801fbd8f88SHong Zhang   ierr = MPI_Comm_split(comm,color,subrank,&subcomm);CHKERRQ(ierr);
811fbd8f88SHong Zhang 
821fbd8f88SHong Zhang   j = 0; duprank = 0;
831fbd8f88SHong Zhang   for (i=0; i<nsubcomm; i++){
841fbd8f88SHong Zhang     if (j == color){
851fbd8f88SHong Zhang       duprank += subrank;
861fbd8f88SHong Zhang       break;
871fbd8f88SHong Zhang     }
881fbd8f88SHong Zhang     duprank += subsize[i]; j++;
891fbd8f88SHong Zhang   }
901fbd8f88SHong Zhang 
911fbd8f88SHong Zhang   /* create dupcomm with same size as comm, but its rank, duprank, maps subcomm's contiguously into dupcomm */
921fbd8f88SHong Zhang   ierr = MPI_Comm_split(comm,0,duprank,&dupcomm);CHKERRQ(ierr);
931fbd8f88SHong Zhang   ierr = PetscFree(subsize);CHKERRQ(ierr);
941fbd8f88SHong Zhang 
951fbd8f88SHong Zhang   ierr = PetscNew(PetscSubcomm,&psubcomm_tmp);CHKERRQ(ierr);
961fbd8f88SHong Zhang   psubcomm_tmp->parent    = comm;
971fbd8f88SHong Zhang   psubcomm_tmp->dupparent = dupcomm;
981fbd8f88SHong Zhang   psubcomm_tmp->comm      = subcomm;
991fbd8f88SHong Zhang   psubcomm_tmp->n         = nsubcomm;
1001fbd8f88SHong Zhang   psubcomm_tmp->color     = color;
1011fbd8f88SHong Zhang   *psubcomm = psubcomm_tmp;
1021fbd8f88SHong Zhang   PetscFunctionReturn(0);
1031fbd8f88SHong Zhang }
1041fbd8f88SHong Zhang 
1051fbd8f88SHong Zhang typedef struct {
1064b9ad928SBarry Smith   PC           pc;                   /* actual preconditioner used on each processor */
1073f457be1SHong Zhang   Vec          xsub,ysub;            /* vectors of a subcommunicator to hold parallel vectors of pc->comm */
1083f457be1SHong Zhang   Vec          xdup,ydup;            /* parallel vector that congregates xsub or ysub facilitating vector scattering */
109b3804887SHong Zhang   Mat          pmats;                /* matrix and optional preconditioner matrix belong to a subcommunicator */
1103f457be1SHong Zhang   VecScatter   scatterin,scatterout; /* scatter used to move all values to each processor group (subcommunicator) */
1114b9ad928SBarry Smith   PetscTruth   useparallelmat;
1121fbd8f88SHong Zhang   PetscSubcomm *psubcomm;
1131fbd8f88SHong Zhang   PetscInt     nsubcomm;           /* num of data structure PetscSubcomm */
1144b9ad928SBarry Smith } PC_Redundant;
1154b9ad928SBarry Smith 
1164b9ad928SBarry Smith #undef __FUNCT__
1174b9ad928SBarry Smith #define __FUNCT__ "PCView_Redundant"
1186849ba73SBarry Smith static PetscErrorCode PCView_Redundant(PC pc,PetscViewer viewer)
1194b9ad928SBarry Smith {
1204b9ad928SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
121dfbe8321SBarry Smith   PetscErrorCode ierr;
12213f74950SBarry Smith   PetscMPIInt    rank;
12332077d6dSBarry Smith   PetscTruth     iascii,isstring;
124a47c9f9aSHong Zhang   PetscViewer    sviewer,subviewer;
1251fbd8f88SHong Zhang   PetscInt       color = red->psubcomm->color;
1264b9ad928SBarry Smith 
1274b9ad928SBarry Smith   PetscFunctionBegin;
1284b9ad928SBarry Smith   ierr = MPI_Comm_rank(pc->comm,&rank);CHKERRQ(ierr);
12932077d6dSBarry Smith   ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr);
1304b9ad928SBarry Smith   ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);CHKERRQ(ierr);
13132077d6dSBarry Smith   if (iascii) {
132a98ce0f4SHong Zhang     ierr = PetscViewerASCIIPrintf(viewer,"  Redundant solver preconditioner: First PC (color=0) follows\n");CHKERRQ(ierr);
133a98ce0f4SHong Zhang     ierr = PetscViewerGetSubcomm(viewer,red->pc->comm,&subviewer);CHKERRQ(ierr);
134a98ce0f4SHong Zhang     if (!color) { /* only view first redundant pc */
1354b9ad928SBarry Smith       ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
136a47c9f9aSHong Zhang       ierr = PCView(red->pc,subviewer);CHKERRQ(ierr);
1374b9ad928SBarry Smith       ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
1384b9ad928SBarry Smith     }
139a98ce0f4SHong Zhang     ierr = PetscViewerRestoreSubcomm(viewer,red->pc->comm,&subviewer);CHKERRQ(ierr);
140a98ce0f4SHong Zhang   } else if (isstring) { /* not test it yet! */
1414b9ad928SBarry Smith     ierr = PetscViewerStringSPrintf(viewer," Redundant solver preconditioner");CHKERRQ(ierr);
1424b9ad928SBarry Smith     ierr = PetscViewerGetSingleton(viewer,&sviewer);CHKERRQ(ierr);
1434b9ad928SBarry Smith     if (!rank) {
1444b9ad928SBarry Smith       ierr = PCView(red->pc,sviewer);CHKERRQ(ierr);
1454b9ad928SBarry Smith     }
1464b9ad928SBarry Smith     ierr = PetscViewerRestoreSingleton(viewer,&sviewer);CHKERRQ(ierr);
1474b9ad928SBarry Smith   } else {
14879a5c55eSBarry Smith     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PC redundant",((PetscObject)viewer)->type_name);
1494b9ad928SBarry Smith   }
1504b9ad928SBarry Smith   PetscFunctionReturn(0);
1514b9ad928SBarry Smith }
1524b9ad928SBarry Smith 
153b9147fbbSdalcinl #include "include/private/matimpl.h"        /*I "petscmat.h" I*/
1544b9ad928SBarry Smith #undef __FUNCT__
1554b9ad928SBarry Smith #define __FUNCT__ "PCSetUp_Redundant"
1566849ba73SBarry Smith static PetscErrorCode PCSetUp_Redundant(PC pc)
1574b9ad928SBarry Smith {
1584b9ad928SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
159dfbe8321SBarry Smith   PetscErrorCode ierr;
1603f457be1SHong Zhang   PetscInt       mstart,mend,mlocal,m;
16113f74950SBarry Smith   PetscMPIInt    size;
1624b9ad928SBarry Smith   MatReuse       reuse = MAT_INITIAL_MATRIX;
1634b9ad928SBarry Smith   MatStructure   str   = DIFFERENT_NONZERO_PATTERN;
1641fbd8f88SHong Zhang   MPI_Comm       comm = pc->comm;
16523ce1328SBarry Smith   Vec            vec;
1663f457be1SHong Zhang   PetscInt       mlocal_sub;
1673f457be1SHong Zhang   PetscMPIInt    subsize,subrank;
1681fbd8f88SHong Zhang   PetscInt       rstart_sub,rend_sub,mloc_sub,nsubcomm;
1691fbd8f88SHong Zhang   const char     *prefix;
1703f457be1SHong Zhang 
1714b9ad928SBarry Smith   PetscFunctionBegin;
17223ce1328SBarry Smith   ierr = MatGetVecs(pc->pmat,&vec,0);CHKERRQ(ierr);
17323ce1328SBarry Smith   ierr = VecGetSize(vec,&m);CHKERRQ(ierr);
1741fbd8f88SHong Zhang 
1754b9ad928SBarry Smith   if (!pc->setupcalled) {
1761fbd8f88SHong Zhang     ierr = PetscSubcommCreate(comm,red->nsubcomm,&red->psubcomm);CHKERRQ(ierr);
1771fbd8f88SHong Zhang     ierr = PetscLogObjectMemory(pc,sizeof(PetscSubcomm));CHKERRQ(ierr);
1781fbd8f88SHong Zhang 
1791fbd8f88SHong Zhang     /* create a new PC that processors in each subcomm have copy of */
1801fbd8f88SHong Zhang     MPI_Comm subcomm = red->psubcomm->comm;
1811fbd8f88SHong Zhang     ierr = PCCreate(subcomm,&red->pc);CHKERRQ(ierr);
1821fbd8f88SHong Zhang     ierr = PCSetType(red->pc,PCLU);CHKERRQ(ierr);
1831fbd8f88SHong Zhang     ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr);
1841fbd8f88SHong Zhang     ierr = PCSetOptionsPrefix(red->pc,prefix);CHKERRQ(ierr);
1851fbd8f88SHong Zhang     ierr = PCAppendOptionsPrefix(red->pc,"redundant_");CHKERRQ(ierr);
1861fbd8f88SHong Zhang     ierr = PCSetFromOptions(red->pc);CHKERRQ(ierr);
1871fbd8f88SHong Zhang 
1883f457be1SHong Zhang     /* create working vectors xsub/ysub and xdup/ydup */
18923ce1328SBarry Smith     ierr = VecGetLocalSize(vec,&mlocal);CHKERRQ(ierr);
1903f457be1SHong Zhang     ierr = VecGetOwnershipRange(vec,&mstart,&mend);CHKERRQ(ierr);
1914b9ad928SBarry Smith 
1923f457be1SHong Zhang     /* get local size of xsub/ysub */
1931fbd8f88SHong Zhang     ierr = MPI_Comm_size(subcomm,&subsize);CHKERRQ(ierr);
1941fbd8f88SHong Zhang     ierr = MPI_Comm_rank(subcomm,&subrank);CHKERRQ(ierr);
1951fbd8f88SHong Zhang     rstart_sub = pc->pmat->rmap.range[red->psubcomm->n*subrank]; /* rstart in xsub/ysub */
1963f457be1SHong Zhang     if (subrank+1 < subsize){
1971fbd8f88SHong Zhang       rend_sub = pc->pmat->rmap.range[red->psubcomm->n*(subrank+1)];
1983f457be1SHong Zhang     } else {
1993f457be1SHong Zhang       rend_sub = m;
2003f457be1SHong Zhang     }
2013f457be1SHong Zhang     mloc_sub = rend_sub - rstart_sub;
2021fbd8f88SHong Zhang     ierr = VecCreateMPI(subcomm,mloc_sub,PETSC_DECIDE,&red->ysub);CHKERRQ(ierr);
2033f457be1SHong Zhang     /* create xsub with empty local arrays, because xdup's arrays will be placed into it */
2041fbd8f88SHong Zhang     ierr = VecCreateMPIWithArray(subcomm,mloc_sub,PETSC_DECIDE,PETSC_NULL,&red->xsub);CHKERRQ(ierr);
2053f457be1SHong Zhang 
2063f457be1SHong Zhang     /* create xdup and ydup. ydup has empty local arrays because ysub's arrays will be place into it.
2073f457be1SHong Zhang        Note: we use communicator dupcomm, not pc->comm! */
2081fbd8f88SHong Zhang     ierr = VecCreateMPI(red->psubcomm->dupparent,mloc_sub,PETSC_DECIDE,&red->xdup);CHKERRQ(ierr);
2091fbd8f88SHong Zhang     ierr = VecCreateMPIWithArray(red->psubcomm->dupparent,mloc_sub,PETSC_DECIDE,PETSC_NULL,&red->ydup);CHKERRQ(ierr);
2103f457be1SHong Zhang 
2113f457be1SHong Zhang     /* create vec scatters */
2123f457be1SHong Zhang     if (!red->scatterin){
2133f457be1SHong Zhang       IS       is1,is2;
2143f457be1SHong Zhang       PetscInt *idx1,*idx2,i,j,k;
2151fbd8f88SHong Zhang       ierr = PetscMalloc(2*red->psubcomm->n*mlocal*sizeof(PetscInt),&idx1);CHKERRQ(ierr);
2161fbd8f88SHong Zhang       idx2 = idx1 + red->psubcomm->n*mlocal;
2173f457be1SHong Zhang       j = 0;
2181fbd8f88SHong Zhang       for (k=0; k<red->psubcomm->n; k++){
2193f457be1SHong Zhang         for (i=mstart; i<mend; i++){
2203f457be1SHong Zhang           idx1[j]   = i;
2213f457be1SHong Zhang           idx2[j++] = i + m*k;
2223f457be1SHong Zhang         }
2233f457be1SHong Zhang       }
2241fbd8f88SHong Zhang       ierr = ISCreateGeneral(comm,red->psubcomm->n*mlocal,idx1,&is1);CHKERRQ(ierr);
2251fbd8f88SHong Zhang       ierr = ISCreateGeneral(comm,red->psubcomm->n*mlocal,idx2,&is2);CHKERRQ(ierr);
2263f457be1SHong Zhang       ierr = VecScatterCreate(vec,is1,red->xdup,is2,&red->scatterin);CHKERRQ(ierr);
2273f457be1SHong Zhang       ierr = ISDestroy(is1);CHKERRQ(ierr);
2283f457be1SHong Zhang       ierr = ISDestroy(is2);CHKERRQ(ierr);
2293f457be1SHong Zhang 
2301fbd8f88SHong Zhang       ierr = ISCreateStride(comm,mlocal,mstart+ red->psubcomm->color*m,1,&is1);CHKERRQ(ierr);
2313f457be1SHong Zhang       ierr = ISCreateStride(comm,mlocal,mstart,1,&is2);CHKERRQ(ierr);
2323f457be1SHong Zhang       ierr = VecScatterCreate(red->xdup,is1,vec,is2,&red->scatterout);CHKERRQ(ierr);
2333f457be1SHong Zhang       ierr = ISDestroy(is1);CHKERRQ(ierr);
2343f457be1SHong Zhang       ierr = ISDestroy(is2);CHKERRQ(ierr);
2353f457be1SHong Zhang       ierr = PetscFree(idx1);CHKERRQ(ierr);
2364b9ad928SBarry Smith     }
2374b9ad928SBarry Smith   }
23823ce1328SBarry Smith   ierr = VecDestroy(vec);CHKERRQ(ierr);
2394b9ad928SBarry Smith 
2404b9ad928SBarry Smith   /* if pmatrix set by user is sequential then we do not need to gather the parallel matrix */
2413f457be1SHong Zhang   ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
2424b9ad928SBarry Smith   if (size == 1) {
2434b9ad928SBarry Smith     red->useparallelmat = PETSC_FALSE;
2444b9ad928SBarry Smith   }
2454b9ad928SBarry Smith 
2464b9ad928SBarry Smith   if (red->useparallelmat) {
2474b9ad928SBarry Smith     if (pc->setupcalled == 1 && pc->flag == DIFFERENT_NONZERO_PATTERN) {
2484b9ad928SBarry Smith       /* destroy old matrices */
2494b9ad928SBarry Smith       if (red->pmats) {
250b3804887SHong Zhang         ierr = MatDestroy(red->pmats);CHKERRQ(ierr);
2514b9ad928SBarry Smith       }
2524b9ad928SBarry Smith     } else if (pc->setupcalled == 1) {
2534b9ad928SBarry Smith       reuse = MAT_REUSE_MATRIX;
2544b9ad928SBarry Smith       str   = SAME_NONZERO_PATTERN;
2554b9ad928SBarry Smith     }
2564b9ad928SBarry Smith 
2573f457be1SHong Zhang     /* grab the parallel matrix and put it into processors of a subcomminicator */
258f664ae05SHong Zhang     /*--------------------------------------------------------------------------*/
259f664ae05SHong Zhang     ierr = VecGetLocalSize(red->ysub,&mlocal_sub);CHKERRQ(ierr);
260*69db28dcSHong Zhang     ierr = MatGetRedundantMatrix(pc->pmat,red->psubcomm->n,red->psubcomm->comm,mlocal_sub,reuse,&red->pmats);CHKERRQ(ierr);
261*69db28dcSHong Zhang 
2623f457be1SHong Zhang     /* tell PC of the subcommunicator its operators */
263b3804887SHong Zhang     ierr = PCSetOperators(red->pc,red->pmats,red->pmats,str);CHKERRQ(ierr);
2644b9ad928SBarry Smith   } else {
2654b9ad928SBarry Smith     ierr = PCSetOperators(red->pc,pc->mat,pc->pmat,pc->flag);CHKERRQ(ierr);
2664b9ad928SBarry Smith   }
2674b9ad928SBarry Smith   ierr = PCSetFromOptions(red->pc);CHKERRQ(ierr);
2684b9ad928SBarry Smith   ierr = PCSetUp(red->pc);CHKERRQ(ierr);
2694b9ad928SBarry Smith   PetscFunctionReturn(0);
2704b9ad928SBarry Smith }
2714b9ad928SBarry Smith 
2724b9ad928SBarry Smith #undef __FUNCT__
2734b9ad928SBarry Smith #define __FUNCT__ "PCApply_Redundant"
2746849ba73SBarry Smith static PetscErrorCode PCApply_Redundant(PC pc,Vec x,Vec y)
2754b9ad928SBarry Smith {
2764b9ad928SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
277dfbe8321SBarry Smith   PetscErrorCode ierr;
2783f457be1SHong Zhang   PetscScalar    *array;
2794b9ad928SBarry Smith 
2804b9ad928SBarry Smith   PetscFunctionBegin;
2813f457be1SHong Zhang   /* scatter x to xdup */
2823f457be1SHong Zhang   ierr = VecScatterBegin(x,red->xdup,INSERT_VALUES,SCATTER_FORWARD,red->scatterin);CHKERRQ(ierr);
2833f457be1SHong Zhang   ierr = VecScatterEnd(x,red->xdup,INSERT_VALUES,SCATTER_FORWARD,red->scatterin);CHKERRQ(ierr);
2843f457be1SHong Zhang 
2853f457be1SHong Zhang   /* place xdup's local array into xsub */
2863f457be1SHong Zhang   ierr = VecGetArray(red->xdup,&array);CHKERRQ(ierr);
2873f457be1SHong Zhang   ierr = VecPlaceArray(red->xsub,(const PetscScalar*)array);CHKERRQ(ierr);
2884b9ad928SBarry Smith 
2894b9ad928SBarry Smith   /* apply preconditioner on each processor */
2903f457be1SHong Zhang   ierr = PCApply(red->pc,red->xsub,red->ysub);CHKERRQ(ierr);
2913f457be1SHong Zhang   ierr = VecResetArray(red->xsub);CHKERRQ(ierr);
2923f457be1SHong Zhang   ierr = VecRestoreArray(red->xdup,&array);CHKERRQ(ierr);
2934b9ad928SBarry Smith 
2943f457be1SHong Zhang   /* place ysub's local array into ydup */
2953f457be1SHong Zhang   ierr = VecGetArray(red->ysub,&array);CHKERRQ(ierr);
2963f457be1SHong Zhang   ierr = VecPlaceArray(red->ydup,(const PetscScalar*)array);CHKERRQ(ierr);
2973f457be1SHong Zhang 
2983f457be1SHong Zhang   /* scatter ydup to y */
2993f457be1SHong Zhang   ierr = VecScatterBegin(red->ydup,y,INSERT_VALUES,SCATTER_FORWARD,red->scatterout);CHKERRQ(ierr);
3003f457be1SHong Zhang   ierr = VecScatterEnd(red->ydup,y,INSERT_VALUES,SCATTER_FORWARD,red->scatterout);CHKERRQ(ierr);
3013f457be1SHong Zhang   ierr = VecResetArray(red->ydup);CHKERRQ(ierr);
3023f457be1SHong Zhang   ierr = VecRestoreArray(red->ysub,&array);CHKERRQ(ierr);
3034b9ad928SBarry Smith   PetscFunctionReturn(0);
3044b9ad928SBarry Smith }
3054b9ad928SBarry Smith 
3064b9ad928SBarry Smith #undef __FUNCT__
3074b9ad928SBarry Smith #define __FUNCT__ "PCDestroy_Redundant"
3086849ba73SBarry Smith static PetscErrorCode PCDestroy_Redundant(PC pc)
3094b9ad928SBarry Smith {
3104b9ad928SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
311dfbe8321SBarry Smith   PetscErrorCode ierr;
3124b9ad928SBarry Smith 
3134b9ad928SBarry Smith   PetscFunctionBegin;
3144b9ad928SBarry Smith   if (red->scatterin)  {ierr = VecScatterDestroy(red->scatterin);CHKERRQ(ierr);}
3154b9ad928SBarry Smith   if (red->scatterout) {ierr = VecScatterDestroy(red->scatterout);CHKERRQ(ierr);}
3163f457be1SHong Zhang   if (red->ysub)       {ierr = VecDestroy(red->ysub);CHKERRQ(ierr);}
3173f457be1SHong Zhang   if (red->xsub)       {ierr = VecDestroy(red->xsub);CHKERRQ(ierr);}
3183f457be1SHong Zhang   if (red->xdup)       {ierr = VecDestroy(red->xdup);CHKERRQ(ierr);}
3193f457be1SHong Zhang   if (red->ydup)       {ierr = VecDestroy(red->ydup);CHKERRQ(ierr);}
320b3804887SHong Zhang   if (red->pmats) {
321b3804887SHong Zhang     ierr = MatDestroy(red->pmats);CHKERRQ(ierr);
3223f457be1SHong Zhang   }
3231fbd8f88SHong Zhang   ierr = PetscSubcommDestroy(red->psubcomm);CHKERRQ(ierr);
3244b9ad928SBarry Smith   ierr = PCDestroy(red->pc);CHKERRQ(ierr);
3254b9ad928SBarry Smith   ierr = PetscFree(red);CHKERRQ(ierr);
3264b9ad928SBarry Smith   PetscFunctionReturn(0);
3274b9ad928SBarry Smith }
3284b9ad928SBarry Smith 
3294b9ad928SBarry Smith #undef __FUNCT__
3304b9ad928SBarry Smith #define __FUNCT__ "PCSetFromOptions_Redundant"
3316849ba73SBarry Smith static PetscErrorCode PCSetFromOptions_Redundant(PC pc)
3324b9ad928SBarry Smith {
333a98ce0f4SHong Zhang   PetscErrorCode ierr;
334a98ce0f4SHong Zhang   PC_Redundant   *red = (PC_Redundant*)pc->data;
3351fbd8f88SHong Zhang   PetscMPIInt    size;
336a98ce0f4SHong Zhang 
3374b9ad928SBarry Smith   PetscFunctionBegin;
338a98ce0f4SHong Zhang   ierr = PetscOptionsHead("Redundant options");CHKERRQ(ierr);
339*69db28dcSHong Zhang   ierr = PetscOptionsInt("-pc_redundant_number_comm","Number of subcommunicators","PCRedundantSetNumComm",red->nsubcomm,&red->nsubcomm,0);CHKERRQ(ierr);
340a98ce0f4SHong Zhang   ierr = PetscOptionsTail();CHKERRQ(ierr);
3414b9ad928SBarry Smith   PetscFunctionReturn(0);
3424b9ad928SBarry Smith }
3434b9ad928SBarry Smith 
3444b9ad928SBarry Smith EXTERN_C_BEGIN
3454b9ad928SBarry Smith #undef __FUNCT__
3464b9ad928SBarry Smith #define __FUNCT__ "PCRedundantSetScatter_Redundant"
347dba47a55SKris Buschelman PetscErrorCode PETSCKSP_DLLEXPORT PCRedundantSetScatter_Redundant(PC pc,VecScatter in,VecScatter out)
3484b9ad928SBarry Smith {
3494b9ad928SBarry Smith   PC_Redundant   *red = (PC_Redundant*)pc->data;
350dfbe8321SBarry Smith   PetscErrorCode ierr;
3514b9ad928SBarry Smith 
3524b9ad928SBarry Smith   PetscFunctionBegin;
3534b9ad928SBarry Smith   red->scatterin  = in;
3544b9ad928SBarry Smith   red->scatterout = out;
3554b9ad928SBarry Smith   ierr = PetscObjectReference((PetscObject)in);CHKERRQ(ierr);
3564b9ad928SBarry Smith   ierr = PetscObjectReference((PetscObject)out);CHKERRQ(ierr);
3574b9ad928SBarry Smith   PetscFunctionReturn(0);
3584b9ad928SBarry Smith }
3594b9ad928SBarry Smith EXTERN_C_END
3604b9ad928SBarry Smith 
3614b9ad928SBarry Smith #undef __FUNCT__
3624b9ad928SBarry Smith #define __FUNCT__ "PCRedundantSetScatter"
3634b9ad928SBarry Smith /*@
3644b9ad928SBarry Smith    PCRedundantSetScatter - Sets the scatter used to copy values into the
3654b9ad928SBarry Smith      redundant local solve and the scatter to move them back into the global
3664b9ad928SBarry Smith      vector.
3674b9ad928SBarry Smith 
3684b9ad928SBarry Smith    Collective on PC
3694b9ad928SBarry Smith 
3704b9ad928SBarry Smith    Input Parameters:
3714b9ad928SBarry Smith +  pc - the preconditioner context
3724b9ad928SBarry Smith .  in - the scatter to move the values in
3734b9ad928SBarry Smith -  out - the scatter to move them out
3744b9ad928SBarry Smith 
3754b9ad928SBarry Smith    Level: advanced
3764b9ad928SBarry Smith 
3774b9ad928SBarry Smith .keywords: PC, redundant solve
3784b9ad928SBarry Smith @*/
379dba47a55SKris Buschelman PetscErrorCode PETSCKSP_DLLEXPORT PCRedundantSetScatter(PC pc,VecScatter in,VecScatter out)
3804b9ad928SBarry Smith {
381dfbe8321SBarry Smith   PetscErrorCode ierr,(*f)(PC,VecScatter,VecScatter);
3824b9ad928SBarry Smith 
3834b9ad928SBarry Smith   PetscFunctionBegin;
3844482741eSBarry Smith   PetscValidHeaderSpecific(pc,PC_COOKIE,1);
3854482741eSBarry Smith   PetscValidHeaderSpecific(in,VEC_SCATTER_COOKIE,2);
3864482741eSBarry Smith   PetscValidHeaderSpecific(out,VEC_SCATTER_COOKIE,3);
3874b9ad928SBarry Smith   ierr = PetscObjectQueryFunction((PetscObject)pc,"PCRedundantSetScatter_C",(void (**)(void))&f);CHKERRQ(ierr);
3884b9ad928SBarry Smith   if (f) {
3894b9ad928SBarry Smith     ierr = (*f)(pc,in,out);CHKERRQ(ierr);
3904b9ad928SBarry Smith   }
3914b9ad928SBarry Smith   PetscFunctionReturn(0);
3924b9ad928SBarry Smith }
3934b9ad928SBarry Smith 
3944b9ad928SBarry Smith EXTERN_C_BEGIN
3954b9ad928SBarry Smith #undef __FUNCT__
3964b9ad928SBarry Smith #define __FUNCT__ "PCRedundantGetPC_Redundant"
397dba47a55SKris Buschelman PetscErrorCode PETSCKSP_DLLEXPORT PCRedundantGetPC_Redundant(PC pc,PC *innerpc)
3984b9ad928SBarry Smith {
3994b9ad928SBarry Smith   PC_Redundant *red = (PC_Redundant*)pc->data;
4004b9ad928SBarry Smith 
4014b9ad928SBarry Smith   PetscFunctionBegin;
4024b9ad928SBarry Smith   *innerpc = red->pc;
4034b9ad928SBarry Smith   PetscFunctionReturn(0);
4044b9ad928SBarry Smith }
4054b9ad928SBarry Smith EXTERN_C_END
4064b9ad928SBarry Smith 
4074b9ad928SBarry Smith #undef __FUNCT__
4084b9ad928SBarry Smith #define __FUNCT__ "PCRedundantGetPC"
4094b9ad928SBarry Smith /*@
4104b9ad928SBarry Smith    PCRedundantGetPC - Gets the sequential PC created by the redundant PC.
4114b9ad928SBarry Smith 
4124b9ad928SBarry Smith    Not Collective
4134b9ad928SBarry Smith 
4144b9ad928SBarry Smith    Input Parameter:
4154b9ad928SBarry Smith .  pc - the preconditioner context
4164b9ad928SBarry Smith 
4174b9ad928SBarry Smith    Output Parameter:
4184b9ad928SBarry Smith .  innerpc - the sequential PC
4194b9ad928SBarry Smith 
4204b9ad928SBarry Smith    Level: advanced
4214b9ad928SBarry Smith 
4224b9ad928SBarry Smith .keywords: PC, redundant solve
4234b9ad928SBarry Smith @*/
424dba47a55SKris Buschelman PetscErrorCode PETSCKSP_DLLEXPORT PCRedundantGetPC(PC pc,PC *innerpc)
4254b9ad928SBarry Smith {
426dfbe8321SBarry Smith   PetscErrorCode ierr,(*f)(PC,PC*);
4274b9ad928SBarry Smith 
4284b9ad928SBarry Smith   PetscFunctionBegin;
4294482741eSBarry Smith   PetscValidHeaderSpecific(pc,PC_COOKIE,1);
4304482741eSBarry Smith   PetscValidPointer(innerpc,2);
4314b9ad928SBarry Smith   ierr = PetscObjectQueryFunction((PetscObject)pc,"PCRedundantGetPC_C",(void (**)(void))&f);CHKERRQ(ierr);
4324b9ad928SBarry Smith   if (f) {
4334b9ad928SBarry Smith     ierr = (*f)(pc,innerpc);CHKERRQ(ierr);
4344b9ad928SBarry Smith   }
4354b9ad928SBarry Smith   PetscFunctionReturn(0);
4364b9ad928SBarry Smith }
4374b9ad928SBarry Smith 
4384b9ad928SBarry Smith EXTERN_C_BEGIN
4394b9ad928SBarry Smith #undef __FUNCT__
4404b9ad928SBarry Smith #define __FUNCT__ "PCRedundantGetOperators_Redundant"
441dba47a55SKris Buschelman PetscErrorCode PETSCKSP_DLLEXPORT PCRedundantGetOperators_Redundant(PC pc,Mat *mat,Mat *pmat)
4424b9ad928SBarry Smith {
4434b9ad928SBarry Smith   PC_Redundant *red = (PC_Redundant*)pc->data;
4444b9ad928SBarry Smith 
4454b9ad928SBarry Smith   PetscFunctionBegin;
446b3804887SHong Zhang   if (mat)  *mat  = red->pmats;
447b3804887SHong Zhang   if (pmat) *pmat = red->pmats;
4484b9ad928SBarry Smith   PetscFunctionReturn(0);
4494b9ad928SBarry Smith }
4504b9ad928SBarry Smith EXTERN_C_END
4514b9ad928SBarry Smith 
4524b9ad928SBarry Smith #undef __FUNCT__
4534b9ad928SBarry Smith #define __FUNCT__ "PCRedundantGetOperators"
4544b9ad928SBarry Smith /*@
4554b9ad928SBarry Smith    PCRedundantGetOperators - gets the sequential matrix and preconditioner matrix
4564b9ad928SBarry Smith 
4574b9ad928SBarry Smith    Not Collective
4584b9ad928SBarry Smith 
4594b9ad928SBarry Smith    Input Parameter:
4604b9ad928SBarry Smith .  pc - the preconditioner context
4614b9ad928SBarry Smith 
4624b9ad928SBarry Smith    Output Parameters:
4634b9ad928SBarry Smith +  mat - the matrix
4644b9ad928SBarry Smith -  pmat - the (possibly different) preconditioner matrix
4654b9ad928SBarry Smith 
4664b9ad928SBarry Smith    Level: advanced
4674b9ad928SBarry Smith 
4684b9ad928SBarry Smith .keywords: PC, redundant solve
4694b9ad928SBarry Smith @*/
470dba47a55SKris Buschelman PetscErrorCode PETSCKSP_DLLEXPORT PCRedundantGetOperators(PC pc,Mat *mat,Mat *pmat)
4714b9ad928SBarry Smith {
472dfbe8321SBarry Smith   PetscErrorCode ierr,(*f)(PC,Mat*,Mat*);
4734b9ad928SBarry Smith 
4744b9ad928SBarry Smith   PetscFunctionBegin;
4754482741eSBarry Smith   PetscValidHeaderSpecific(pc,PC_COOKIE,1);
4764482741eSBarry Smith   if (mat)  PetscValidPointer(mat,2);
4774482741eSBarry Smith   if (pmat) PetscValidPointer(pmat,3);
4784b9ad928SBarry Smith   ierr = PetscObjectQueryFunction((PetscObject)pc,"PCRedundantGetOperators_C",(void (**)(void))&f);CHKERRQ(ierr);
4794b9ad928SBarry Smith   if (f) {
4804b9ad928SBarry Smith     ierr = (*f)(pc,mat,pmat);CHKERRQ(ierr);
4814b9ad928SBarry Smith   }
4824b9ad928SBarry Smith   PetscFunctionReturn(0);
4834b9ad928SBarry Smith }
4844b9ad928SBarry Smith 
4854b9ad928SBarry Smith /* -------------------------------------------------------------------------------------*/
48637a17b4dSBarry Smith /*MC
487a98ce0f4SHong Zhang      PCREDUNDANT - Runs a preconditioner for the entire problem on subgroups of processors
48837a17b4dSBarry Smith 
48937a17b4dSBarry Smith      Options for the redundant preconditioners can be set with -redundant_pc_xxx
49037a17b4dSBarry Smith 
49109391456SBarry Smith   Options Database:
49209391456SBarry Smith .  -pc_redundant_number_comm - number of sub communicators to use
49309391456SBarry Smith 
49437a17b4dSBarry Smith    Level: intermediate
49537a17b4dSBarry Smith 
49637a17b4dSBarry Smith .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PCRedundantSetScatter(),
49737a17b4dSBarry Smith            PCRedundantGetPC(), PCRedundantGetOperators()
49837a17b4dSBarry Smith M*/
49937a17b4dSBarry Smith 
5004b9ad928SBarry Smith EXTERN_C_BEGIN
5014b9ad928SBarry Smith #undef __FUNCT__
5024b9ad928SBarry Smith #define __FUNCT__ "PCCreate_Redundant"
503dba47a55SKris Buschelman PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_Redundant(PC pc)
5044b9ad928SBarry Smith {
505dfbe8321SBarry Smith   PetscErrorCode ierr;
5064b9ad928SBarry Smith   PC_Redundant   *red;
507*69db28dcSHong Zhang   PetscMPIInt    size;
5083f457be1SHong Zhang 
5094b9ad928SBarry Smith   PetscFunctionBegin;
5104b9ad928SBarry Smith   ierr = PetscNew(PC_Redundant,&red);CHKERRQ(ierr);
51152e6d16bSBarry Smith   ierr = PetscLogObjectMemory(pc,sizeof(PC_Redundant));CHKERRQ(ierr);
512*69db28dcSHong Zhang   ierr = MPI_Comm_size(pc->comm,&size);CHKERRQ(ierr);
513*69db28dcSHong Zhang   red->nsubcomm       = size;
5144b9ad928SBarry Smith   red->useparallelmat = PETSC_TRUE;
5151fbd8f88SHong Zhang   pc->data            = (void*)red;
5164b9ad928SBarry Smith 
5174b9ad928SBarry Smith   pc->ops->apply           = PCApply_Redundant;
5184b9ad928SBarry Smith   pc->ops->applytranspose  = 0;
5194b9ad928SBarry Smith   pc->ops->setup           = PCSetUp_Redundant;
5204b9ad928SBarry Smith   pc->ops->destroy         = PCDestroy_Redundant;
5214b9ad928SBarry Smith   pc->ops->setfromoptions  = PCSetFromOptions_Redundant;
5224b9ad928SBarry Smith   pc->ops->view            = PCView_Redundant;
5234b9ad928SBarry Smith   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCRedundantSetScatter_C","PCRedundantSetScatter_Redundant",
5244b9ad928SBarry Smith                     PCRedundantSetScatter_Redundant);CHKERRQ(ierr);
5254b9ad928SBarry Smith   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCRedundantGetPC_C","PCRedundantGetPC_Redundant",
5264b9ad928SBarry Smith                     PCRedundantGetPC_Redundant);CHKERRQ(ierr);
5274b9ad928SBarry Smith   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCRedundantGetOperators_C","PCRedundantGetOperators_Redundant",
5284b9ad928SBarry Smith                     PCRedundantGetOperators_Redundant);CHKERRQ(ierr);
5294b9ad928SBarry Smith   PetscFunctionReturn(0);
5304b9ad928SBarry Smith }
5314b9ad928SBarry Smith EXTERN_C_END
532