xref: /petsc/src/ksp/pc/impls/redundant/redundant.c (revision ddc54837807eb96c5bc9ef34be37eed03061befe)
1 
2 /*
3   This file defines a "solve the problem redundantly on each subgroup of processor" preconditioner.
4 */
5 #include <petsc-private/pcimpl.h>
6 #include <petscksp.h>           /*I "petscksp.h" I*/
7 
8 typedef struct {
9   KSP          ksp;
10   PC           pc;                   /* actual preconditioner used on each processor */
11   Vec          xsub,ysub;            /* vectors of a subcommunicator to hold parallel vectors of PetscObjectComm((PetscObject)pc) */
12   Vec          xdup,ydup;            /* parallel vector that congregates xsub or ysub facilitating vector scattering */
13   Mat          pmats;                /* matrix and optional preconditioner matrix belong to a subcommunicator */
14   VecScatter   scatterin,scatterout; /* scatter used to move all values to each processor group (subcommunicator) */
15   PetscBool    useparallelmat;
16   PetscSubcomm psubcomm;
17   PetscInt     nsubcomm;           /* num of data structure PetscSubcomm */
18 } PC_Redundant;
19 
20 #undef __FUNCT__
21 #define __FUNCT__ "PCView_Redundant"
22 static PetscErrorCode PCView_Redundant(PC pc,PetscViewer viewer)
23 {
24   PC_Redundant   *red = (PC_Redundant*)pc->data;
25   PetscErrorCode ierr;
26   PetscBool      iascii,isstring;
27   PetscViewer    subviewer;
28 
29   PetscFunctionBegin;
30   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
31   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);CHKERRQ(ierr);
32   if (iascii) {
33     if (!red->psubcomm) {
34       ierr = PetscViewerASCIIPrintf(viewer,"  Redundant preconditioner: Not yet setup\n");CHKERRQ(ierr);
35     } else {
36       ierr = PetscViewerASCIIPrintf(viewer,"  Redundant preconditioner: First (color=0) of %D PCs follows\n",red->nsubcomm);CHKERRQ(ierr);
37       ierr = PetscViewerGetSubcomm(viewer,((PetscObject)red->pc)->comm,&subviewer);CHKERRQ(ierr);
38       if (!red->psubcomm->color) { /* only view first redundant pc */
39         ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
40         ierr = KSPView(red->ksp,subviewer);CHKERRQ(ierr);
41         ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
42       }
43       ierr = PetscViewerRestoreSubcomm(viewer,((PetscObject)red->pc)->comm,&subviewer);CHKERRQ(ierr);
44     }
45   } else if (isstring) {
46     ierr = PetscViewerStringSPrintf(viewer," Redundant solver preconditioner");CHKERRQ(ierr);
47   }
48   PetscFunctionReturn(0);
49 }
50 
51 extern PetscErrorCode MatGetRedundantMatrix_MPIAIJ_interlaced(Mat,PetscInt,PetscSubcomm,MatReuse,Mat*); /* rm later! */
52 
53 #undef __FUNCT__
54 #define __FUNCT__ "PCSetUp_Redundant"
55 static PetscErrorCode PCSetUp_Redundant(PC pc)
56 {
57   PC_Redundant   *red = (PC_Redundant*)pc->data;
58   PetscErrorCode ierr;
59   PetscInt       mstart,mend,mlocal,M,mlocal_sub,rstart_sub,rend_sub,mloc_sub;
60   PetscMPIInt    size;
61   MatReuse       reuse = MAT_INITIAL_MATRIX;
62   MatStructure   str   = DIFFERENT_NONZERO_PATTERN;
63   MPI_Comm       comm,subcomm;
64   Vec            x;
65   PetscMPIInt    subsize,subrank;
66   const char     *prefix;
67   const PetscInt *range;
68 
69   PetscFunctionBegin;
70   ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
71   ierr = MatGetVecs(pc->pmat,&x,0);CHKERRQ(ierr);
72   ierr = VecGetSize(x,&M);CHKERRQ(ierr);
73 
74   /* if pmatrix set by user is sequential then we do not need to gather the parallel matrix */
75   ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
76   if (size == 1) red->useparallelmat = PETSC_FALSE;
77 
78   if (!pc->setupcalled) {
79     if (!red->psubcomm) {
80       ierr = PetscSubcommCreate(comm,&red->psubcomm);CHKERRQ(ierr);
81       ierr = PetscSubcommSetNumber(red->psubcomm,red->nsubcomm);CHKERRQ(ierr);
82       ierr = PetscSubcommSetType(red->psubcomm,PETSC_SUBCOMM_INTERLACED);CHKERRQ(ierr);
83       ierr = PetscLogObjectMemory(pc,sizeof(PetscSubcomm));CHKERRQ(ierr);
84 
85       /* create a new PC that processors in each subcomm have copy of */
86       subcomm = red->psubcomm->comm;
87 
88       ierr = KSPCreate(subcomm,&red->ksp);CHKERRQ(ierr);
89       ierr = PetscObjectIncrementTabLevel((PetscObject)red->ksp,(PetscObject)pc,1);CHKERRQ(ierr);
90       ierr = PetscLogObjectParent(pc,red->ksp);CHKERRQ(ierr);
91       ierr = KSPSetType(red->ksp,KSPPREONLY);CHKERRQ(ierr);
92       ierr = KSPGetPC(red->ksp,&red->pc);CHKERRQ(ierr);
93       ierr = PCSetType(red->pc,PCLU);CHKERRQ(ierr);
94 
95       ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr);
96       ierr = KSPSetOptionsPrefix(red->ksp,prefix);CHKERRQ(ierr);
97       ierr = KSPAppendOptionsPrefix(red->ksp,"redundant_");CHKERRQ(ierr);
98     } else subcomm = red->psubcomm->comm;
99 
100     /* create working vectors xsub/ysub and xdup/ydup */
101     ierr = VecGetLocalSize(x,&mlocal);CHKERRQ(ierr);
102     ierr = VecGetOwnershipRange(x,&mstart,&mend);CHKERRQ(ierr);
103 
104     /* get local size of xsub/ysub */
105     ierr = MPI_Comm_size(subcomm,&subsize);CHKERRQ(ierr);
106     ierr = MPI_Comm_rank(subcomm,&subrank);CHKERRQ(ierr);
107     ierr = MatGetOwnershipRanges(pc->pmat,&range);CHKERRQ(ierr);
108     rstart_sub = range[red->psubcomm->n*subrank]; /* rstart in xsub/ysub */
109     if (subrank+1 < subsize) rend_sub = range[red->psubcomm->n*(subrank+1)];
110     else rend_sub = M;
111 
112     mloc_sub = rend_sub - rstart_sub;
113     ierr     = VecCreateMPI(subcomm,mloc_sub,PETSC_DECIDE,&red->ysub);CHKERRQ(ierr);
114     /* create xsub with empty local arrays, because xdup's arrays will be placed into it */
115     ierr = VecCreateMPIWithArray(subcomm,1,mloc_sub,PETSC_DECIDE,NULL,&red->xsub);CHKERRQ(ierr);
116 
117     /* create xdup and ydup. ydup has empty local arrays because ysub's arrays will be place into it.
118        Note: we use communicator dupcomm, not PetscObjectComm((PetscObject)pc)! */
119     ierr = VecCreateMPI(red->psubcomm->dupparent,mloc_sub,PETSC_DECIDE,&red->xdup);CHKERRQ(ierr);
120     ierr = VecCreateMPIWithArray(red->psubcomm->dupparent,1,mloc_sub,PETSC_DECIDE,NULL,&red->ydup);CHKERRQ(ierr);
121 
122     /* create vec scatters */
123     if (!red->scatterin) {
124       IS       is1,is2;
125       PetscInt *idx1,*idx2,i,j,k;
126 
127       ierr = PetscMalloc2(red->psubcomm->n*mlocal,PetscInt,&idx1,red->psubcomm->n*mlocal,PetscInt,&idx2);CHKERRQ(ierr);
128       j    = 0;
129       for (k=0; k<red->psubcomm->n; k++) {
130         for (i=mstart; i<mend; i++) {
131           idx1[j]   = i;
132           idx2[j++] = i + M*k;
133         }
134       }
135       ierr = ISCreateGeneral(comm,red->psubcomm->n*mlocal,idx1,PETSC_COPY_VALUES,&is1);CHKERRQ(ierr);
136       ierr = ISCreateGeneral(comm,red->psubcomm->n*mlocal,idx2,PETSC_COPY_VALUES,&is2);CHKERRQ(ierr);
137       ierr = VecScatterCreate(x,is1,red->xdup,is2,&red->scatterin);CHKERRQ(ierr);
138       ierr = ISDestroy(&is1);CHKERRQ(ierr);
139       ierr = ISDestroy(&is2);CHKERRQ(ierr);
140 
141       ierr = ISCreateStride(comm,mlocal,mstart+ red->psubcomm->color*M,1,&is1);CHKERRQ(ierr);
142       ierr = ISCreateStride(comm,mlocal,mstart,1,&is2);CHKERRQ(ierr);
143       ierr = VecScatterCreate(red->xdup,is1,x,is2,&red->scatterout);CHKERRQ(ierr);
144       ierr = ISDestroy(&is1);CHKERRQ(ierr);
145       ierr = ISDestroy(&is2);CHKERRQ(ierr);
146       ierr = PetscFree2(idx1,idx2);CHKERRQ(ierr);
147     }
148   }
149   ierr = VecDestroy(&x);CHKERRQ(ierr);
150 
151   if (red->useparallelmat) {
152     if (pc->setupcalled == 1 && pc->flag == DIFFERENT_NONZERO_PATTERN) {
153       /* destroy old matrices */
154       ierr = MatDestroy(&red->pmats);CHKERRQ(ierr);
155     } else if (pc->setupcalled == 1) {
156       reuse = MAT_REUSE_MATRIX;
157       str   = SAME_NONZERO_PATTERN;
158     }
159 
160     /* grab the parallel matrix and put it into processors of a subcomminicator */
161     /*--------------------------------------------------------------------------*/
162     ierr = VecGetLocalSize(red->ysub,&mlocal_sub);CHKERRQ(ierr);
163     ierr = MatGetRedundantMatrix_MPIAIJ_interlaced(pc->pmat,red->psubcomm->n,red->psubcomm,reuse,&red->pmats);CHKERRQ(ierr);
164     /* tell PC of the subcommunicator its operators */
165     ierr = KSPSetOperators(red->ksp,red->pmats,red->pmats,str);CHKERRQ(ierr);
166   } else {
167     ierr = KSPSetOperators(red->ksp,pc->mat,pc->pmat,pc->flag);CHKERRQ(ierr);
168   }
169   if (pc->setfromoptionscalled) {
170     ierr = KSPSetFromOptions(red->ksp);CHKERRQ(ierr);
171   }
172   ierr = KSPSetUp(red->ksp);CHKERRQ(ierr);
173   PetscFunctionReturn(0);
174 }
175 
176 #undef __FUNCT__
177 #define __FUNCT__ "PCApply_Redundant"
178 static PetscErrorCode PCApply_Redundant(PC pc,Vec x,Vec y)
179 {
180   PC_Redundant   *red = (PC_Redundant*)pc->data;
181   PetscErrorCode ierr;
182   PetscScalar    *array;
183 
184   PetscFunctionBegin;
185   if (!red->useparallelmat) {
186     ierr = KSPSolve(red->ksp,x,y);CHKERRQ(ierr);
187     PetscFunctionReturn(0);
188   }
189 
190   /* scatter x to xdup */
191   ierr = VecScatterBegin(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
192   ierr = VecScatterEnd(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
193 
194   /* place xdup's local array into xsub */
195   ierr = VecGetArray(red->xdup,&array);CHKERRQ(ierr);
196   ierr = VecPlaceArray(red->xsub,(const PetscScalar*)array);CHKERRQ(ierr);
197 
198   /* apply preconditioner on each processor */
199   ierr = KSPSolve(red->ksp,red->xsub,red->ysub);CHKERRQ(ierr);
200   ierr = VecResetArray(red->xsub);CHKERRQ(ierr);
201   ierr = VecRestoreArray(red->xdup,&array);CHKERRQ(ierr);
202 
203   /* place ysub's local array into ydup */
204   ierr = VecGetArray(red->ysub,&array);CHKERRQ(ierr);
205   ierr = VecPlaceArray(red->ydup,(const PetscScalar*)array);CHKERRQ(ierr);
206 
207   /* scatter ydup to y */
208   ierr = VecScatterBegin(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
209   ierr = VecScatterEnd(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
210   ierr = VecResetArray(red->ydup);CHKERRQ(ierr);
211   ierr = VecRestoreArray(red->ysub,&array);CHKERRQ(ierr);
212   PetscFunctionReturn(0);
213 }
214 
215 #undef __FUNCT__
216 #define __FUNCT__ "PCReset_Redundant"
217 static PetscErrorCode PCReset_Redundant(PC pc)
218 {
219   PC_Redundant   *red = (PC_Redundant*)pc->data;
220   PetscErrorCode ierr;
221 
222   PetscFunctionBegin;
223   ierr = VecScatterDestroy(&red->scatterin);CHKERRQ(ierr);
224   ierr = VecScatterDestroy(&red->scatterout);CHKERRQ(ierr);
225   ierr = VecDestroy(&red->ysub);CHKERRQ(ierr);
226   ierr = VecDestroy(&red->xsub);CHKERRQ(ierr);
227   ierr = VecDestroy(&red->xdup);CHKERRQ(ierr);
228   ierr = VecDestroy(&red->ydup);CHKERRQ(ierr);
229   ierr = MatDestroy(&red->pmats);CHKERRQ(ierr);
230   if (red->ksp) {ierr = KSPReset(red->ksp);CHKERRQ(ierr);}
231   PetscFunctionReturn(0);
232 }
233 
234 #undef __FUNCT__
235 #define __FUNCT__ "PCDestroy_Redundant"
236 static PetscErrorCode PCDestroy_Redundant(PC pc)
237 {
238   PC_Redundant   *red = (PC_Redundant*)pc->data;
239   PetscErrorCode ierr;
240 
241   PetscFunctionBegin;
242   ierr = PCReset_Redundant(pc);CHKERRQ(ierr);
243   ierr = KSPDestroy(&red->ksp);CHKERRQ(ierr);
244   ierr = PetscSubcommDestroy(&red->psubcomm);CHKERRQ(ierr);
245   ierr = PetscFree(pc->data);CHKERRQ(ierr);
246   PetscFunctionReturn(0);
247 }
248 
249 #undef __FUNCT__
250 #define __FUNCT__ "PCSetFromOptions_Redundant"
251 static PetscErrorCode PCSetFromOptions_Redundant(PC pc)
252 {
253   PetscErrorCode ierr;
254   PC_Redundant   *red = (PC_Redundant*)pc->data;
255 
256   PetscFunctionBegin;
257   ierr = PetscOptionsHead("Redundant options");CHKERRQ(ierr);
258   ierr = PetscOptionsInt("-pc_redundant_number","Number of redundant pc","PCRedundantSetNumber",red->nsubcomm,&red->nsubcomm,0);CHKERRQ(ierr);
259   ierr = PetscOptionsTail();CHKERRQ(ierr);
260   PetscFunctionReturn(0);
261 }
262 
263 #undef __FUNCT__
264 #define __FUNCT__ "PCRedundantSetNumber_Redundant"
265 static PetscErrorCode  PCRedundantSetNumber_Redundant(PC pc,PetscInt nreds)
266 {
267   PC_Redundant *red = (PC_Redundant*)pc->data;
268 
269   PetscFunctionBegin;
270   red->nsubcomm = nreds;
271   PetscFunctionReturn(0);
272 }
273 
274 #undef __FUNCT__
275 #define __FUNCT__ "PCRedundantSetNumber"
276 /*@
277    PCRedundantSetNumber - Sets the number of redundant preconditioner contexts.
278 
279    Logically Collective on PC
280 
281    Input Parameters:
282 +  pc - the preconditioner context
283 -  nredundant - number of redundant preconditioner contexts; for example if you are using 64 MPI processes and
284                               use an nredundant of 4 there will be 4 parallel solves each on 16 = 64/4 processes.
285 
286    Level: advanced
287 
288 .keywords: PC, redundant solve
289 @*/
290 PetscErrorCode  PCRedundantSetNumber(PC pc,PetscInt nredundant)
291 {
292   PetscErrorCode ierr;
293 
294   PetscFunctionBegin;
295   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
296   if (nredundant <= 0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "num of redundant pc %D must be positive",nredundant);
297   ierr = PetscTryMethod(pc,"PCRedundantSetNumber_C",(PC,PetscInt),(pc,nredundant));CHKERRQ(ierr);
298   PetscFunctionReturn(0);
299 }
300 
301 #undef __FUNCT__
302 #define __FUNCT__ "PCRedundantSetScatter_Redundant"
303 static PetscErrorCode  PCRedundantSetScatter_Redundant(PC pc,VecScatter in,VecScatter out)
304 {
305   PC_Redundant   *red = (PC_Redundant*)pc->data;
306   PetscErrorCode ierr;
307 
308   PetscFunctionBegin;
309   ierr = PetscObjectReference((PetscObject)in);CHKERRQ(ierr);
310   ierr = VecScatterDestroy(&red->scatterin);CHKERRQ(ierr);
311 
312   red->scatterin  = in;
313 
314   ierr = PetscObjectReference((PetscObject)out);CHKERRQ(ierr);
315   ierr = VecScatterDestroy(&red->scatterout);CHKERRQ(ierr);
316   red->scatterout = out;
317   PetscFunctionReturn(0);
318 }
319 
320 #undef __FUNCT__
321 #define __FUNCT__ "PCRedundantSetScatter"
322 /*@
323    PCRedundantSetScatter - Sets the scatter used to copy values into the
324      redundant local solve and the scatter to move them back into the global
325      vector.
326 
327    Logically Collective on PC
328 
329    Input Parameters:
330 +  pc - the preconditioner context
331 .  in - the scatter to move the values in
332 -  out - the scatter to move them out
333 
334    Level: advanced
335 
336 .keywords: PC, redundant solve
337 @*/
338 PetscErrorCode  PCRedundantSetScatter(PC pc,VecScatter in,VecScatter out)
339 {
340   PetscErrorCode ierr;
341 
342   PetscFunctionBegin;
343   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
344   PetscValidHeaderSpecific(in,VEC_SCATTER_CLASSID,2);
345   PetscValidHeaderSpecific(out,VEC_SCATTER_CLASSID,3);
346   ierr = PetscTryMethod(pc,"PCRedundantSetScatter_C",(PC,VecScatter,VecScatter),(pc,in,out));CHKERRQ(ierr);
347   PetscFunctionReturn(0);
348 }
349 
350 #undef __FUNCT__
351 #define __FUNCT__ "PCRedundantGetKSP_Redundant"
352 static PetscErrorCode  PCRedundantGetKSP_Redundant(PC pc,KSP *innerksp)
353 {
354   PetscErrorCode ierr;
355   PC_Redundant   *red = (PC_Redundant*)pc->data;
356   MPI_Comm       comm,subcomm;
357   const char     *prefix;
358 
359   PetscFunctionBegin;
360   if (!red->psubcomm) {
361     ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
362     ierr = PetscSubcommCreate(comm,&red->psubcomm);CHKERRQ(ierr);
363     ierr = PetscSubcommSetNumber(red->psubcomm,red->nsubcomm);CHKERRQ(ierr);
364     ierr = PetscSubcommSetType(red->psubcomm,PETSC_SUBCOMM_INTERLACED);CHKERRQ(ierr);
365     ierr = PetscLogObjectMemory(pc,sizeof(PetscSubcomm));CHKERRQ(ierr);
366 
367     /* create a new PC that processors in each subcomm have copy of */
368     subcomm = red->psubcomm->comm;
369 
370     ierr = KSPCreate(subcomm,&red->ksp);CHKERRQ(ierr);
371     ierr = PetscObjectIncrementTabLevel((PetscObject)red->ksp,(PetscObject)pc,1);CHKERRQ(ierr);
372     ierr = PetscLogObjectParent(pc,red->ksp);CHKERRQ(ierr);
373     ierr = KSPSetType(red->ksp,KSPPREONLY);CHKERRQ(ierr);
374     ierr = KSPGetPC(red->ksp,&red->pc);CHKERRQ(ierr);
375     ierr = PCSetType(red->pc,PCLU);CHKERRQ(ierr);
376 
377     ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr);
378     ierr = KSPSetOptionsPrefix(red->ksp,prefix);CHKERRQ(ierr);
379     ierr = KSPAppendOptionsPrefix(red->ksp,"redundant_");CHKERRQ(ierr);
380   }
381   *innerksp = red->ksp;
382   PetscFunctionReturn(0);
383 }
384 
385 #undef __FUNCT__
386 #define __FUNCT__ "PCRedundantGetKSP"
387 /*@
388    PCRedundantGetKSP - Gets the less parallel KSP created by the redundant PC.
389 
390    Not Collective
391 
392    Input Parameter:
393 .  pc - the preconditioner context
394 
395    Output Parameter:
396 .  innerksp - the KSP on the smaller set of processes
397 
398    Level: advanced
399 
400 .keywords: PC, redundant solve
401 @*/
402 PetscErrorCode  PCRedundantGetKSP(PC pc,KSP *innerksp)
403 {
404   PetscErrorCode ierr;
405 
406   PetscFunctionBegin;
407   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
408   PetscValidPointer(innerksp,2);
409   ierr = PetscTryMethod(pc,"PCRedundantGetKSP_C",(PC,KSP*),(pc,innerksp));CHKERRQ(ierr);
410   PetscFunctionReturn(0);
411 }
412 
413 #undef __FUNCT__
414 #define __FUNCT__ "PCRedundantGetOperators_Redundant"
415 static PetscErrorCode  PCRedundantGetOperators_Redundant(PC pc,Mat *mat,Mat *pmat)
416 {
417   PC_Redundant *red = (PC_Redundant*)pc->data;
418 
419   PetscFunctionBegin;
420   if (mat)  *mat  = red->pmats;
421   if (pmat) *pmat = red->pmats;
422   PetscFunctionReturn(0);
423 }
424 
425 #undef __FUNCT__
426 #define __FUNCT__ "PCRedundantGetOperators"
427 /*@
428    PCRedundantGetOperators - gets the sequential matrix and preconditioner matrix
429 
430    Not Collective
431 
432    Input Parameter:
433 .  pc - the preconditioner context
434 
435    Output Parameters:
436 +  mat - the matrix
437 -  pmat - the (possibly different) preconditioner matrix
438 
439    Level: advanced
440 
441 .keywords: PC, redundant solve
442 @*/
443 PetscErrorCode  PCRedundantGetOperators(PC pc,Mat *mat,Mat *pmat)
444 {
445   PetscErrorCode ierr;
446 
447   PetscFunctionBegin;
448   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
449   if (mat)  PetscValidPointer(mat,2);
450   if (pmat) PetscValidPointer(pmat,3);
451   ierr = PetscTryMethod(pc,"PCRedundantGetOperators_C",(PC,Mat*,Mat*),(pc,mat,pmat));CHKERRQ(ierr);
452   PetscFunctionReturn(0);
453 }
454 
455 /* -------------------------------------------------------------------------------------*/
456 /*MC
457      PCREDUNDANT - Runs a KSP solver with preconditioner for the entire problem on subgroups of processors
458 
459      Options for the redundant preconditioners can be set with -redundant_pc_xxx for the redundant KSP with -redundant_ksp_xxx
460 
461   Options Database:
462 .  -pc_redundant_number <n> - number of redundant solves, for example if you are using 64 MPI processes and
463                               use an n of 4 there will be 4 parallel solves each on 16 = 64/4 processes.
464 
465    Level: intermediate
466 
467    Notes: The default KSP is preonly and the default PC is LU.
468 
469    Developer Notes: Note that PCSetInitialGuessNonzero()  is not used by this class but likely should be.
470 
471 .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PCRedundantSetScatter(),
472            PCRedundantGetKSP(), PCRedundantGetOperators(), PCRedundantSetNumber()
473 M*/
474 
475 #undef __FUNCT__
476 #define __FUNCT__ "PCCreate_Redundant"
477 PETSC_EXTERN PetscErrorCode PCCreate_Redundant(PC pc)
478 {
479   PetscErrorCode ierr;
480   PC_Redundant   *red;
481   PetscMPIInt    size;
482 
483   PetscFunctionBegin;
484   ierr = PetscNewLog(pc,PC_Redundant,&red);CHKERRQ(ierr);
485   ierr = MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);CHKERRQ(ierr);
486 
487   red->nsubcomm       = size;
488   red->useparallelmat = PETSC_TRUE;
489   pc->data            = (void*)red;
490 
491   pc->ops->apply          = PCApply_Redundant;
492   pc->ops->applytranspose = 0;
493   pc->ops->setup          = PCSetUp_Redundant;
494   pc->ops->destroy        = PCDestroy_Redundant;
495   pc->ops->reset          = PCReset_Redundant;
496   pc->ops->setfromoptions = PCSetFromOptions_Redundant;
497   pc->ops->view           = PCView_Redundant;
498 
499   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCRedundantSetScatter_C",PCRedundantSetScatter_Redundant);CHKERRQ(ierr);
500   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCRedundantSetNumber_C",PCRedundantSetNumber_Redundant);CHKERRQ(ierr);
501   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCRedundantGetKSP_C",PCRedundantGetKSP_Redundant);CHKERRQ(ierr);
502   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCRedundantGetOperators_C",PCRedundantGetOperators_Redundant);CHKERRQ(ierr);
503   PetscFunctionReturn(0);
504 }
505 
506