xref: /petsc/src/ksp/pc/impls/fieldsplit/fieldsplit.c (revision 5e45c6bbc0badcad264181c3e3af51a3b424ec19)
1 #define PETSCKSP_DLL
2 
3 /*
4 
5 */
6 #include "private/pcimpl.h"     /*I "petscpc.h" I*/
7 
8 typedef struct _PC_FieldSplitLink *PC_FieldSplitLink;
9 struct _PC_FieldSplitLink {
10   KSP               ksp;
11   Vec               x,y;
12   PetscInt          nfields;
13   PetscInt          *fields;
14   VecScatter        sctx;
15   PC_FieldSplitLink next,previous;
16 };
17 
18 typedef struct {
19   PCCompositeType   type;              /* additive or multiplicative */
20   PetscTruth        defaultsplit;
21   PetscInt          bs;
22   PetscInt          nsplits,*csize;
23   Vec               *x,*y,w1,w2;
24   Mat               *pmat;
25   IS                *is,*cis;
26   PC_FieldSplitLink head;
27 } PC_FieldSplit;
28 
29 #undef __FUNCT__
30 #define __FUNCT__ "PCView_FieldSplit"
31 static PetscErrorCode PCView_FieldSplit(PC pc,PetscViewer viewer)
32 {
33   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
34   PetscErrorCode    ierr;
35   PetscTruth        iascii;
36   PetscInt          i,j;
37   PC_FieldSplitLink ilink = jac->head;
38 
39   PetscFunctionBegin;
40   ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr);
41   if (iascii) {
42     ierr = PetscViewerASCIIPrintf(viewer,"  FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs);CHKERRQ(ierr);
43     ierr = PetscViewerASCIIPrintf(viewer,"  Solver info for each split is in the following KSP objects:\n");CHKERRQ(ierr);
44     ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
45     for (i=0; i<jac->nsplits; i++) {
46       ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);CHKERRQ(ierr);
47       ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr);
48       for (j=0; j<ilink->nfields; j++) {
49         if (j > 0) {
50           ierr = PetscViewerASCIIPrintf(viewer,",");CHKERRQ(ierr);
51         }
52         ierr = PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);CHKERRQ(ierr);
53       }
54       ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr);
55       ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr);
56       ierr = KSPView(ilink->ksp,viewer);CHKERRQ(ierr);
57       ilink = ilink->next;
58     }
59     ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
60   } else {
61     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCFieldSplit",((PetscObject)viewer)->type_name);
62   }
63   PetscFunctionReturn(0);
64 }
65 
66 #undef __FUNCT__
67 #define __FUNCT__ "PCFieldSplitSetDefaults"
68 static PetscErrorCode PCFieldSplitSetDefaults(PC pc)
69 {
70   PC_FieldSplit     *jac  = (PC_FieldSplit*)pc->data;
71   PetscErrorCode    ierr;
72   PC_FieldSplitLink ilink = jac->head;
73   PetscInt          i;
74   PetscTruth        flg = PETSC_FALSE,*fields;
75 
76   PetscFunctionBegin;
77   ierr = PetscOptionsGetTruth(pc->prefix,"-pc_fieldsplit_default",&flg,PETSC_NULL);CHKERRQ(ierr);
78   if (!ilink || flg) {
79     ierr = PetscInfo(pc,"Using default splitting of fields\n");CHKERRQ(ierr);
80     if (jac->bs <= 0) {
81       ierr   = MatGetBlockSize(pc->pmat,&jac->bs);CHKERRQ(ierr);
82     }
83     ierr = PetscMalloc(jac->bs*sizeof(PetscTruth),&fields);CHKERRQ(ierr);
84     ierr = PetscMemzero(fields,jac->bs*sizeof(PetscTruth));CHKERRQ(ierr);
85     while (ilink) {
86       for (i=0; i<ilink->nfields; i++) {
87         fields[ilink->fields[i]] = PETSC_TRUE;
88       }
89       ilink = ilink->next;
90     }
91     jac->defaultsplit = PETSC_TRUE;
92     for (i=0; i<jac->bs; i++) {
93       if (!fields[i]) {
94 	ierr = PCFieldSplitSetFields(pc,1,&i);CHKERRQ(ierr);
95       } else {
96         jac->defaultsplit = PETSC_FALSE;
97       }
98     }
99   }
100   PetscFunctionReturn(0);
101 }
102 
103 
104 #undef __FUNCT__
105 #define __FUNCT__ "PCSetUp_FieldSplit"
106 static PetscErrorCode PCSetUp_FieldSplit(PC pc)
107 {
108   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
109   PetscErrorCode    ierr;
110   PC_FieldSplitLink ilink;
111   PetscInt          i,nsplit,ccsize;
112   MatStructure      flag = pc->flag;
113 
114   PetscFunctionBegin;
115   ierr   = PCFieldSplitSetDefaults(pc);CHKERRQ(ierr);
116   nsplit = jac->nsplits;
117   ilink  = jac->head;
118 
119   /* get the matrices for each split */
120   if (!jac->is) {
121     PetscInt rstart,rend,nslots,bs;
122 
123     bs     = jac->bs;
124     ierr   = MatGetOwnershipRange(pc->pmat,&rstart,&rend);CHKERRQ(ierr);
125     ierr   = MatGetLocalSize(pc->pmat,PETSC_NULL,&ccsize);CHKERRQ(ierr);
126     nslots = (rend - rstart)/bs;
127     ierr   = PetscMalloc(nsplit*sizeof(IS),&jac->is);CHKERRQ(ierr);
128     ierr   = PetscMalloc(nsplit*sizeof(IS),&jac->cis);CHKERRQ(ierr);
129     ierr   = PetscMalloc(nsplit*sizeof(PetscInt),&jac->csize);CHKERRQ(ierr);
130     for (i=0; i<nsplit; i++) {
131       if (jac->defaultsplit) {
132 	ierr     = ISCreateStride(pc->comm,nslots,rstart+i,nsplit,&jac->is[i]);CHKERRQ(ierr);
133         jac->csize[i] = ccsize/nsplit;
134       } else {
135         PetscInt   *ii,j,k,nfields = ilink->nfields,*fields = ilink->fields;
136         PetscTruth sorted;
137         ierr = PetscMalloc(ilink->nfields*nslots*sizeof(PetscInt),&ii);CHKERRQ(ierr);
138         for (j=0; j<nslots; j++) {
139           for (k=0; k<nfields; k++) {
140             ii[nfields*j + k] = rstart + bs*j + fields[k];
141           }
142         }
143 	ierr = ISCreateGeneral(pc->comm,nslots*nfields,ii,&jac->is[i]);CHKERRQ(ierr);
144         jac->csize[i] = (ccsize/bs)*ilink->nfields;
145         ierr = ISSorted(jac->is[i],&sorted);CHKERRQ(ierr);
146         if (!sorted) SETERRQ(PETSC_ERR_USER,"Fields must be sorted when creating split");
147         ierr = PetscFree(ii);CHKERRQ(ierr);
148         ilink = ilink->next;
149       }
150       ierr = ISAllGather(jac->is[i],&jac->cis[i]);CHKERRQ(ierr);
151     }
152   }
153 
154   if (!jac->pmat) {
155     ierr = PetscMalloc(nsplit*sizeof(Mat),&jac->pmat);CHKERRQ(ierr);
156     for (i=0; i<nsplit; i++) {
157       ierr = MatGetSubMatrix(pc->pmat,jac->is[i],jac->cis[i],jac->csize[i],MAT_INITIAL_MATRIX,&jac->pmat[i]);CHKERRQ(ierr);
158     }
159   } else {
160     for (i=0; i<nsplit; i++) {
161       ierr = MatGetSubMatrix(pc->pmat,jac->is[i],jac->cis[i],jac->csize[i],MAT_REUSE_MATRIX,&jac->pmat[i]);CHKERRQ(ierr);
162     }
163   }
164 
165   /* set up the individual PCs */
166   i    = 0;
167   ilink = jac->head;
168   while (ilink) {
169     ierr = KSPSetOperators(ilink->ksp,jac->pmat[i],jac->pmat[i],flag);CHKERRQ(ierr);
170     ierr = KSPSetFromOptions(ilink->ksp);CHKERRQ(ierr);
171     ierr = KSPSetUp(ilink->ksp);CHKERRQ(ierr);
172     i++;
173     ilink = ilink->next;
174   }
175 
176   /* create work vectors for each split */
177   if (!jac->x) {
178     Vec xtmp;
179     ierr = PetscMalloc2(nsplit,Vec,&jac->x,nsplit,Vec,&jac->y);CHKERRQ(ierr);
180     ilink = jac->head;
181     for (i=0; i<nsplit; i++) {
182       Vec *vl,*vr;
183 
184       ierr      = KSPGetVecs(ilink->ksp,1,&vr,1,&vl);CHKERRQ(ierr);
185       ilink->x  = *vr;
186       ilink->y  = *vl;
187       ierr      = PetscFree(vr);CHKERRQ(ierr);
188       ierr      = PetscFree(vl);CHKERRQ(ierr);
189       jac->x[i] = ilink->x;
190       jac->y[i] = ilink->y;
191       ilink     = ilink->next;
192     }
193     /* compute scatter contexts needed by multiplicative versions and non-default splits */
194 
195     ilink = jac->head;
196     ierr = MatGetVecs(pc->pmat,&xtmp,PETSC_NULL);CHKERRQ(ierr);
197     for (i=0; i<nsplit; i++) {
198       ierr = VecScatterCreate(xtmp,jac->is[i],jac->x[i],PETSC_NULL,&ilink->sctx);CHKERRQ(ierr);
199       ilink = ilink->next;
200     }
201     ierr = VecDestroy(xtmp);CHKERRQ(ierr);
202   }
203   PetscFunctionReturn(0);
204 }
205 
206 #define FieldSplitSplitSolveAdd(ilink,xx,yy) \
207     (VecScatterBegin(xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD,ilink->sctx) || \
208      VecScatterEnd(xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD,ilink->sctx) || \
209      KSPSolve(ilink->ksp,ilink->x,ilink->y) || \
210      VecScatterBegin(ilink->y,yy,ADD_VALUES,SCATTER_REVERSE,ilink->sctx) || \
211      VecScatterEnd(ilink->y,yy,ADD_VALUES,SCATTER_REVERSE,ilink->sctx))
212 
213 #undef __FUNCT__
214 #define __FUNCT__ "PCApply_FieldSplit"
215 static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y)
216 {
217   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
218   PetscErrorCode    ierr;
219   PC_FieldSplitLink ilink = jac->head;
220   PetscInt          bs;
221 
222   PetscFunctionBegin;
223   CHKMEMQ;
224   ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr);
225   ierr = VecSetBlockSize(x,jac->bs);CHKERRQ(ierr);
226   ierr = VecSetBlockSize(y,jac->bs);CHKERRQ(ierr);
227 
228   if (jac->type == PC_COMPOSITE_ADDITIVE) {
229     if (jac->defaultsplit) {
230       ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr);
231       while (ilink) {
232 	ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr);
233 	ilink = ilink->next;
234       }
235       ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr);
236     } else {
237       PetscInt    i = 0;
238 
239       ierr = VecSet(y,0.0);CHKERRQ(ierr);
240       while (ilink) {
241         ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr);
242 	ilink = ilink->next;
243 	i++;
244       }
245     }
246   } else {
247     if (!jac->w1) {
248       ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr);
249       ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr);
250     }
251     ierr = VecSet(y,0.0);CHKERRQ(ierr);
252     ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr);
253     while (ilink->next) {
254       ilink = ilink->next;
255       ierr  = MatMult(pc->pmat,y,jac->w1);CHKERRQ(ierr);
256       ierr  = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr);
257       ierr  = FieldSplitSplitSolveAdd(ilink,jac->w2,y);CHKERRQ(ierr);
258     }
259     if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
260       while (ilink->previous) {
261         ilink = ilink->previous;
262         ierr  = MatMult(pc->pmat,y,jac->w1);CHKERRQ(ierr);
263         ierr  = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr);
264         ierr  = FieldSplitSplitSolveAdd(ilink,jac->w2,y);CHKERRQ(ierr);
265       }
266     }
267   }
268   CHKMEMQ;
269   PetscFunctionReturn(0);
270 }
271 
272 #undef __FUNCT__
273 #define __FUNCT__ "PCDestroy_FieldSplit"
274 static PetscErrorCode PCDestroy_FieldSplit(PC pc)
275 {
276   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
277   PetscErrorCode    ierr;
278   PC_FieldSplitLink ilink = jac->head,next;
279 
280   PetscFunctionBegin;
281   while (ilink) {
282     ierr = KSPDestroy(ilink->ksp);CHKERRQ(ierr);
283     if (ilink->x) {ierr = VecDestroy(ilink->x);CHKERRQ(ierr);}
284     if (ilink->y) {ierr = VecDestroy(ilink->y);CHKERRQ(ierr);}
285     if (ilink->sctx) {ierr = VecScatterDestroy(ilink->sctx);CHKERRQ(ierr);}
286     next = ilink->next;
287     ierr = PetscFree2(ilink,ilink->fields);CHKERRQ(ierr);
288     ilink = next;
289   }
290   ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr);
291   if (jac->pmat) {ierr = MatDestroyMatrices(jac->nsplits,&jac->pmat);CHKERRQ(ierr);}
292   if (jac->is) {
293     PetscInt i;
294     for (i=0; i<jac->nsplits; i++) {ierr = ISDestroy(jac->is[i]);CHKERRQ(ierr);}
295     ierr = PetscFree(jac->is);CHKERRQ(ierr);
296   }
297   if (jac->cis) {
298     PetscInt i;
299     for (i=0; i<jac->nsplits; i++) {ierr = ISDestroy(jac->cis[i]);CHKERRQ(ierr);}
300     ierr = PetscFree(jac->cis);CHKERRQ(ierr);
301   }
302   if (jac->w1) {ierr = VecDestroy(jac->w1);CHKERRQ(ierr);}
303   if (jac->w2) {ierr = VecDestroy(jac->w2);CHKERRQ(ierr);}
304   ierr = PetscFree(jac->csize);CHKERRQ(ierr);
305   ierr = PetscFree(jac);CHKERRQ(ierr);
306   PetscFunctionReturn(0);
307 }
308 
309 #undef __FUNCT__
310 #define __FUNCT__ "PCSetFromOptions_FieldSplit"
311 static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc)
312 {
313   PetscErrorCode ierr;
314   PetscInt       i = 0,nfields,*fields,bs;
315   PetscTruth     flg;
316   char           optionname[128];
317   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
318 
319   PetscFunctionBegin;
320   ierr = PetscOptionsHead("FieldSplit options");CHKERRQ(ierr);
321   ierr = PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);CHKERRQ(ierr);
322   if (flg) {
323     ierr = PCFieldSplitSetBlockSize(pc,bs);CHKERRQ(ierr);
324   }
325   if (jac->bs <= 0) {
326     ierr = PCFieldSplitSetBlockSize(pc,1);CHKERRQ(ierr);
327   }
328   ierr = PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&jac->type,&flg);CHKERRQ(ierr);
329   ierr = PetscMalloc(jac->bs*sizeof(PetscInt),&fields);CHKERRQ(ierr);
330   while (PETSC_TRUE) {
331     sprintf(optionname,"-pc_fieldsplit_%d_fields",(int)i++);
332     nfields = jac->bs;
333     ierr    = PetscOptionsIntArray(optionname,"Fields in this split","PCFieldSplitSetFields",fields,&nfields,&flg);CHKERRQ(ierr);
334     if (!flg) break;
335     if (!nfields) SETERRQ(PETSC_ERR_USER,"Cannot list zero fields");
336     ierr = PCFieldSplitSetFields(pc,nfields,fields);CHKERRQ(ierr);
337   }
338   ierr = PetscFree(fields);CHKERRQ(ierr);
339   ierr = PetscOptionsTail();CHKERRQ(ierr);
340   PetscFunctionReturn(0);
341 }
342 
343 /*------------------------------------------------------------------------------------*/
344 
345 EXTERN_C_BEGIN
346 #undef __FUNCT__
347 #define __FUNCT__ "PCFieldSplitSetFields_FieldSplit"
348 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetFields_FieldSplit(PC pc,PetscInt n,PetscInt *fields)
349 {
350   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
351   PetscErrorCode    ierr;
352   PC_FieldSplitLink ilink,next = jac->head;
353   char              prefix[128];
354   PetscInt          i;
355 
356   PetscFunctionBegin;
357   if (n <= 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Negative number of fields requested");
358   for (i=0; i<n; i++) {
359     if (fields[i] >= jac->bs) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs);
360     if (fields[i] < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]);
361   }
362   ierr = PetscMalloc2(1,struct _PC_FieldSplitLink,&ilink,n,PetscInt,&ilink->fields);CHKERRQ(ierr);
363   ierr = PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));CHKERRQ(ierr);
364   ilink->nfields = n;
365   ilink->next    = PETSC_NULL;
366   ierr           = KSPCreate(pc->comm,&ilink->ksp);CHKERRQ(ierr);
367   ierr           = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr);
368 
369   if (pc->prefix) {
370     sprintf(prefix,"%sfieldsplit_%d_",pc->prefix,(int)jac->nsplits);
371   } else {
372     sprintf(prefix,"fieldsplit_%d_",(int)jac->nsplits);
373   }
374   ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr);
375 
376   if (!next) {
377     jac->head       = ilink;
378     ilink->previous = PETSC_NULL;
379   } else {
380     while (next->next) {
381       next = next->next;
382     }
383     next->next      = ilink;
384     ilink->previous = next;
385   }
386   jac->nsplits++;
387   PetscFunctionReturn(0);
388 }
389 EXTERN_C_END
390 
391 
392 EXTERN_C_BEGIN
393 #undef __FUNCT__
394 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit"
395 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp)
396 {
397   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
398   PetscErrorCode    ierr;
399   PetscInt          cnt = 0;
400   PC_FieldSplitLink ilink = jac->head;
401 
402   PetscFunctionBegin;
403   ierr = PetscMalloc(jac->nsplits*sizeof(KSP*),subksp);CHKERRQ(ierr);
404   while (ilink) {
405     (*subksp)[cnt++] = ilink->ksp;
406     ilink = ilink->next;
407   }
408   if (cnt != jac->nsplits) SETERRQ2(PETSC_ERR_PLIB,"Corrupt PCFIELDSPLIT object: number splits in linked list %D in object %D",cnt,jac->nsplits);
409   *n = jac->nsplits;
410   PetscFunctionReturn(0);
411 }
412 EXTERN_C_END
413 
414 #undef __FUNCT__
415 #define __FUNCT__ "PCFieldSplitSetFields"
416 /*@
417     PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner
418 
419     Collective on PC
420 
421     Input Parameters:
422 +   pc  - the preconditioner context
423 .   n - the number of fields in this split
424 .   fields - the fields in this split
425 
426     Level: intermediate
427 
428 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize()
429 
430 @*/
431 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetFields(PC pc,PetscInt n, PetscInt *fields)
432 {
433   PetscErrorCode ierr,(*f)(PC,PetscInt,PetscInt *);
434 
435   PetscFunctionBegin;
436   PetscValidHeaderSpecific(pc,PC_COOKIE,1);
437   ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitSetFields_C",(void (**)(void))&f);CHKERRQ(ierr);
438   if (f) {
439     ierr = (*f)(pc,n,fields);CHKERRQ(ierr);
440   }
441   PetscFunctionReturn(0);
442 }
443 
444 #undef __FUNCT__
445 #define __FUNCT__ "PCFieldSplitSetBlockSize"
446 /*@
447     PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the
448       fieldsplit preconditioner. If not set the matrix block size is used.
449 
450     Collective on PC
451 
452     Input Parameters:
453 +   pc  - the preconditioner context
454 -   bs - the block size
455 
456     Level: intermediate
457 
458 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields()
459 
460 @*/
461 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetBlockSize(PC pc,PetscInt bs)
462 {
463   PetscErrorCode ierr,(*f)(PC,PetscInt);
464 
465   PetscFunctionBegin;
466   PetscValidHeaderSpecific(pc,PC_COOKIE,1);
467   ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",(void (**)(void))&f);CHKERRQ(ierr);
468   if (f) {
469     ierr = (*f)(pc,bs);CHKERRQ(ierr);
470   }
471   PetscFunctionReturn(0);
472 }
473 
474 #undef __FUNCT__
475 #define __FUNCT__ "PCFieldSplitGetSubKSP"
476 /*@C
477    PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits
478 
479    Collective on KSP
480 
481    Input Parameter:
482 .  pc - the preconditioner context
483 
484    Output Parameters:
485 +  n - the number of split
486 -  pc - the array of KSP contexts
487 
488    Note:
489    After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed
490 
491    You must call KSPSetUp() before calling PCFieldSplitGetSubKSP().
492 
493    Level: advanced
494 
495 .seealso: PCFIELDSPLIT
496 @*/
497 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
498 {
499   PetscErrorCode ierr,(*f)(PC,PetscInt*,KSP **);
500 
501   PetscFunctionBegin;
502   PetscValidHeaderSpecific(pc,PC_COOKIE,1);
503   PetscValidIntPointer(n,2);
504   ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",(void (**)(void))&f);CHKERRQ(ierr);
505   if (f) {
506     ierr = (*f)(pc,n,subksp);CHKERRQ(ierr);
507   } else {
508     SETERRQ(PETSC_ERR_ARG_WRONG,"Cannot get subksp for this type of PC");
509   }
510   PetscFunctionReturn(0);
511 }
512 
513 EXTERN_C_BEGIN
514 #undef __FUNCT__
515 #define __FUNCT__ "PCFieldSplitSetType_FieldSplit"
516 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type)
517 {
518   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
519 
520   PetscFunctionBegin;
521   jac->type = type;
522   PetscFunctionReturn(0);
523 }
524 EXTERN_C_END
525 
526 EXTERN_C_BEGIN
527 #undef __FUNCT__
528 #define __FUNCT__ "PCFieldSplitSetBlockSize_FieldSplit"
529 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs)
530 {
531   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
532 
533   PetscFunctionBegin;
534   if (bs < 1) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs);
535   if (jac->bs > 0 && jac->bs != bs) SETERRQ2(PETSC_ERR_ARG_WRONGSTATE,"Cannot change fieldsplit blocksize from %D to %D after it has been set",jac->bs,bs);
536   jac->bs = bs;
537   PetscFunctionReturn(0);
538 }
539 EXTERN_C_END
540 
541 #undef __FUNCT__
542 #define __FUNCT__ "PCFieldSplitSetType"
543 /*@C
544    PCFieldSplitSetType - Sets the type of fieldsplit preconditioner.
545 
546    Collective on PC
547 
548    Input Parameter:
549 .  pc - the preconditioner context
550 .  type - PC_COMPOSITE_ADDITIVE (default), PC_COMPOSITE_MULTIPLICATIVE
551 
552    Options Database Key:
553 .  -pc_fieldsplit_type <type: one of multiplicative, additive, special> - Sets fieldsplit preconditioner type
554 
555    Level: Developer
556 
557 .keywords: PC, set, type, composite preconditioner, additive, multiplicative
558 
559 .seealso: PCCompositeSetType()
560 
561 @*/
562 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetType(PC pc,PCCompositeType type)
563 {
564   PetscErrorCode ierr,(*f)(PC,PCCompositeType);
565 
566   PetscFunctionBegin;
567   PetscValidHeaderSpecific(pc,PC_COOKIE,1);
568   ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitSetType_C",(void (**)(void))&f);CHKERRQ(ierr);
569   if (f) {
570     ierr = (*f)(pc,type);CHKERRQ(ierr);
571   }
572   PetscFunctionReturn(0);
573 }
574 
575 /* -------------------------------------------------------------------------------------*/
576 /*MC
577    PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual
578                   fields or groups of fields
579 
580 
581      To set options on the solvers for each block append -sub_ to all the PC
582         options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1
583 
584      To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP()
585          and set the options directly on the resulting KSP object
586 
587    Level: intermediate
588 
589    Options Database Keys:
590 +   -pc_splitfield_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split
591 .   -pc_splitfield_default - automatically add any fields to additional splits that have not
592                               been supplied explicitly by -pc_splitfield_%d_fields
593 .   -pc_splitfield_block_size <bs> - size of block that defines fields (i.e. there are bs fields)
594 -   -pc_splitfield_type <additive,multiplicative>
595 
596    Concepts: physics based preconditioners
597 
598 .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
599            PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(),PCFieldSplitSetType()
600 M*/
601 
602 EXTERN_C_BEGIN
603 #undef __FUNCT__
604 #define __FUNCT__ "PCCreate_FieldSplit"
605 PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_FieldSplit(PC pc)
606 {
607   PetscErrorCode ierr;
608   PC_FieldSplit  *jac;
609 
610   PetscFunctionBegin;
611   ierr = PetscNew(PC_FieldSplit,&jac);CHKERRQ(ierr);
612   ierr = PetscLogObjectMemory(pc,sizeof(PC_FieldSplit));CHKERRQ(ierr);
613   jac->bs        = -1;
614   jac->nsplits   = 0;
615   jac->type      = PC_COMPOSITE_ADDITIVE;
616 
617   pc->data     = (void*)jac;
618 
619   pc->ops->apply             = PCApply_FieldSplit;
620   pc->ops->setup             = PCSetUp_FieldSplit;
621   pc->ops->destroy           = PCDestroy_FieldSplit;
622   pc->ops->setfromoptions    = PCSetFromOptions_FieldSplit;
623   pc->ops->view              = PCView_FieldSplit;
624   pc->ops->applyrichardson   = 0;
625 
626   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit",
627                     PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr);
628   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","PCFieldSplitSetFields_FieldSplit",
629                     PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr);
630   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetType_C","PCFieldSplitSetType_FieldSplit",
631                     PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr);
632   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetBlockSize_C","PCFieldSplitSetBlockSize_FieldSplit",
633                     PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr);
634   PetscFunctionReturn(0);
635 }
636 EXTERN_C_END
637 
638 
639