Actual source code: fieldsplit.c
petsc-3.9.4 2018-09-11
3: #include <petsc/private/pcimpl.h>
4: #include <petsc/private/kspimpl.h> /* This is needed to provide the appropriate PETSC_EXTERN for KSP_Solve_FS ....*/
5: #include <petscdm.h>
7: const char *const PCFieldSplitSchurPreTypes[] = {"SELF","SELFP","A11","USER","FULL","PCFieldSplitSchurPreType","PC_FIELDSPLIT_SCHUR_PRE_",0};
8: const char *const PCFieldSplitSchurFactTypes[] = {"DIAG","LOWER","UPPER","FULL","PCFieldSplitSchurFactType","PC_FIELDSPLIT_SCHUR_FACT_",0};
10: PetscLogEvent KSP_Solve_FS_0,KSP_Solve_FS_1,KSP_Solve_FS_S,KSP_Solve_FS_U,KSP_Solve_FS_L,KSP_Solve_FS_2,KSP_Solve_FS_3,KSP_Solve_FS_4;
12: typedef struct _PC_FieldSplitLink *PC_FieldSplitLink;
13: struct _PC_FieldSplitLink {
14: KSP ksp;
15: Vec x,y,z;
16: char *splitname;
17: PetscInt nfields;
18: PetscInt *fields,*fields_col;
19: VecScatter sctx;
20: IS is,is_col;
21: PC_FieldSplitLink next,previous;
22: PetscLogEvent event;
23: };
25: typedef struct {
26: PCCompositeType type;
27: PetscBool defaultsplit; /* Flag for a system with a set of 'k' scalar fields with the same layout (and bs = k) */
28: PetscBool splitdefined; /* Flag is set after the splits have been defined, to prevent more splits from being added */
29: PetscInt bs; /* Block size for IS and Mat structures */
30: PetscInt nsplits; /* Number of field divisions defined */
31: Vec *x,*y,w1,w2;
32: Mat *mat; /* The diagonal block for each split */
33: Mat *pmat; /* The preconditioning diagonal block for each split */
34: Mat *Afield; /* The rows of the matrix associated with each split */
35: PetscBool issetup;
37: /* Only used when Schur complement preconditioning is used */
38: Mat B; /* The (0,1) block */
39: Mat C; /* The (1,0) block */
40: Mat schur; /* The Schur complement S = A11 - A10 A00^{-1} A01, the KSP here, kspinner, is H_1 in [El08] */
41: Mat schurp; /* Assembled approximation to S built by MatSchurComplement to be used as a preconditioning matrix when solving with S */
42: Mat schur_user; /* User-provided preconditioning matrix for the Schur complement */
43: PCFieldSplitSchurPreType schurpre; /* Determines which preconditioning matrix is used for the Schur complement */
44: PCFieldSplitSchurFactType schurfactorization;
45: KSP kspschur; /* The solver for S */
46: KSP kspupper; /* The solver for A in the upper diagonal part of the factorization (H_2 in [El08]) */
47: PetscScalar schurscale; /* Scaling factor for the Schur complement solution with DIAG factorization */
49: PC_FieldSplitLink head;
50: PetscBool isrestrict; /* indicates PCFieldSplitRestrictIS() has been last called on this object, hack */
51: PetscBool suboptionsset; /* Indicates that the KSPSetFromOptions() has been called on the sub-KSPs */
52: PetscBool dm_splits; /* Whether to use DMCreateFieldDecomposition() whenever possible */
53: PetscBool diag_use_amat; /* Whether to extract diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */
54: PetscBool offdiag_use_amat; /* Whether to extract off-diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */
55: } PC_FieldSplit;
57: /*
58: Notes: there is no particular reason that pmat, x, and y are stored as arrays in PC_FieldSplit instead of
59: inside PC_FieldSplitLink, just historical. If you want to be able to add new fields after already using the
60: PC you could change this.
61: */
63: /* This helper is so that setting a user-provided preconditioning matrix is orthogonal to choosing to use it. This way the
64: * application-provided FormJacobian can provide this matrix without interfering with the user's (command-line) choices. */
65: static Mat FieldSplitSchurPre(PC_FieldSplit *jac)
66: {
67: switch (jac->schurpre) {
68: case PC_FIELDSPLIT_SCHUR_PRE_SELF: return jac->schur;
69: case PC_FIELDSPLIT_SCHUR_PRE_SELFP: return jac->schurp;
70: case PC_FIELDSPLIT_SCHUR_PRE_A11: return jac->pmat[1];
71: case PC_FIELDSPLIT_SCHUR_PRE_FULL: /* We calculate this and store it in schur_user */
72: case PC_FIELDSPLIT_SCHUR_PRE_USER: /* Use a user-provided matrix if it is given, otherwise diagonal block */
73: default:
74: return jac->schur_user ? jac->schur_user : jac->pmat[1];
75: }
76: }
79: #include <petscdraw.h>
80: static PetscErrorCode PCView_FieldSplit(PC pc,PetscViewer viewer)
81: {
82: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
83: PetscErrorCode ierr;
84: PetscBool iascii,isdraw;
85: PetscInt i,j;
86: PC_FieldSplitLink ilink = jac->head;
89: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
90: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
91: if (iascii) {
92: if (jac->bs > 0) {
93: PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs);
94: } else {
95: PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D\n",PCCompositeTypes[jac->type],jac->nsplits);
96: }
97: if (pc->useAmat) {
98: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");
99: }
100: if (jac->diag_use_amat) {
101: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for diagonal blocks\n");
102: }
103: if (jac->offdiag_use_amat) {
104: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for off-diagonal blocks\n");
105: }
106: PetscViewerASCIIPrintf(viewer," Solver info for each split is in the following KSP objects:\n");
107: PetscViewerASCIIPushTab(viewer);
108: for (i=0; i<jac->nsplits; i++) {
109: if (ilink->fields) {
110: PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);
111: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
112: for (j=0; j<ilink->nfields; j++) {
113: if (j > 0) {
114: PetscViewerASCIIPrintf(viewer,",");
115: }
116: PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);
117: }
118: PetscViewerASCIIPrintf(viewer,"\n");
119: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
120: } else {
121: PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);
122: }
123: KSPView(ilink->ksp,viewer);
124: ilink = ilink->next;
125: }
126: PetscViewerASCIIPopTab(viewer);
127: }
129: if (isdraw) {
130: PetscDraw draw;
131: PetscReal x,y,w,wd;
133: PetscViewerDrawGetDraw(viewer,0,&draw);
134: PetscDrawGetCurrentPoint(draw,&x,&y);
135: w = 2*PetscMin(1.0 - x,x);
136: wd = w/(jac->nsplits + 1);
137: x = x - wd*(jac->nsplits-1)/2.0;
138: for (i=0; i<jac->nsplits; i++) {
139: PetscDrawPushCurrentPoint(draw,x,y);
140: KSPView(ilink->ksp,viewer);
141: PetscDrawPopCurrentPoint(draw);
142: x += wd;
143: ilink = ilink->next;
144: }
145: }
146: return(0);
147: }
149: static PetscErrorCode PCView_FieldSplit_Schur(PC pc,PetscViewer viewer)
150: {
151: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
152: PetscErrorCode ierr;
153: PetscBool iascii,isdraw;
154: PetscInt i,j;
155: PC_FieldSplitLink ilink = jac->head;
156: MatSchurComplementAinvType atype;
159: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
160: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
161: if (iascii) {
162: if (jac->bs > 0) {
163: PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, blocksize = %D, factorization %s\n",jac->bs,PCFieldSplitSchurFactTypes[jac->schurfactorization]);
164: } else {
165: PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, factorization %s\n",PCFieldSplitSchurFactTypes[jac->schurfactorization]);
166: }
167: if (pc->useAmat) {
168: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");
169: }
170: switch (jac->schurpre) {
171: case PC_FIELDSPLIT_SCHUR_PRE_SELF:
172: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from S itself\n");
173: break;
174: case PC_FIELDSPLIT_SCHUR_PRE_SELFP:
175: MatSchurComplementGetAinvType(jac->schur,&atype);
176: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from Sp, an assembled approximation to S, which uses A00's %sdiagonal's inverse\n",atype == MAT_SCHUR_COMPLEMENT_AINV_DIAG ? "" : (atype == MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG ? "block " : "lumped "));break;
177: case PC_FIELDSPLIT_SCHUR_PRE_A11:
178: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n");
179: break;
180: case PC_FIELDSPLIT_SCHUR_PRE_FULL:
181: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from the exact Schur complement\n");
182: break;
183: case PC_FIELDSPLIT_SCHUR_PRE_USER:
184: if (jac->schur_user) {
185: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from user provided matrix\n");
186: } else {
187: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n");
188: }
189: break;
190: default:
191: SETERRQ1(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Invalid Schur preconditioning type: %d", jac->schurpre);
192: }
193: PetscViewerASCIIPrintf(viewer," Split info:\n");
194: PetscViewerASCIIPushTab(viewer);
195: for (i=0; i<jac->nsplits; i++) {
196: if (ilink->fields) {
197: PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);
198: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
199: for (j=0; j<ilink->nfields; j++) {
200: if (j > 0) {
201: PetscViewerASCIIPrintf(viewer,",");
202: }
203: PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);
204: }
205: PetscViewerASCIIPrintf(viewer,"\n");
206: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
207: } else {
208: PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);
209: }
210: ilink = ilink->next;
211: }
212: PetscViewerASCIIPrintf(viewer,"KSP solver for A00 block\n");
213: PetscViewerASCIIPushTab(viewer);
214: if (jac->head) {
215: KSPView(jac->head->ksp,viewer);
216: } else {PetscViewerASCIIPrintf(viewer," not yet available\n");}
217: PetscViewerASCIIPopTab(viewer);
218: if (jac->head && jac->kspupper != jac->head->ksp) {
219: PetscViewerASCIIPrintf(viewer,"KSP solver for upper A00 in upper triangular factor \n");
220: PetscViewerASCIIPushTab(viewer);
221: if (jac->kspupper) {KSPView(jac->kspupper,viewer);}
222: else {PetscViewerASCIIPrintf(viewer," not yet available\n");}
223: PetscViewerASCIIPopTab(viewer);
224: }
225: PetscViewerASCIIPrintf(viewer,"KSP solver for S = A11 - A10 inv(A00) A01 \n");
226: PetscViewerASCIIPushTab(viewer);
227: if (jac->kspschur) {
228: KSPView(jac->kspschur,viewer);
229: } else {
230: PetscViewerASCIIPrintf(viewer," not yet available\n");
231: }
232: PetscViewerASCIIPopTab(viewer);
233: PetscViewerASCIIPopTab(viewer);
234: } else if (isdraw && jac->head) {
235: PetscDraw draw;
236: PetscReal x,y,w,wd,h;
237: PetscInt cnt = 2;
238: char str[32];
240: PetscViewerDrawGetDraw(viewer,0,&draw);
241: PetscDrawGetCurrentPoint(draw,&x,&y);
242: if (jac->kspupper != jac->head->ksp) cnt++;
243: w = 2*PetscMin(1.0 - x,x);
244: wd = w/(cnt + 1);
246: PetscSNPrintf(str,32,"Schur fact. %s",PCFieldSplitSchurFactTypes[jac->schurfactorization]);
247: PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
248: y -= h;
249: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_USER && !jac->schur_user) {
250: PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[PC_FIELDSPLIT_SCHUR_PRE_A11]);
251: } else {
252: PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[jac->schurpre]);
253: }
254: PetscDrawStringBoxed(draw,x+wd*(cnt-1)/2.0,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
255: y -= h;
256: x = x - wd*(cnt-1)/2.0;
258: PetscDrawPushCurrentPoint(draw,x,y);
259: KSPView(jac->head->ksp,viewer);
260: PetscDrawPopCurrentPoint(draw);
261: if (jac->kspupper != jac->head->ksp) {
262: x += wd;
263: PetscDrawPushCurrentPoint(draw,x,y);
264: KSPView(jac->kspupper,viewer);
265: PetscDrawPopCurrentPoint(draw);
266: }
267: x += wd;
268: PetscDrawPushCurrentPoint(draw,x,y);
269: KSPView(jac->kspschur,viewer);
270: PetscDrawPopCurrentPoint(draw);
271: }
272: return(0);
273: }
275: /* Precondition: jac->bs is set to a meaningful value */
276: static PetscErrorCode PCFieldSplitSetRuntimeSplits_Private(PC pc)
277: {
279: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
280: PetscInt i,nfields,*ifields,nfields_col,*ifields_col;
281: PetscBool flg,flg_col;
282: char optionname[128],splitname[8],optionname_col[128];
285: PetscMalloc1(jac->bs,&ifields);
286: PetscMalloc1(jac->bs,&ifields_col);
287: for (i=0,flg=PETSC_TRUE;; i++) {
288: PetscSNPrintf(splitname,sizeof(splitname),"%D",i);
289: PetscSNPrintf(optionname,sizeof(optionname),"-pc_fieldsplit_%D_fields",i);
290: PetscSNPrintf(optionname_col,sizeof(optionname_col),"-pc_fieldsplit_%D_fields_col",i);
291: nfields = jac->bs;
292: nfields_col = jac->bs;
293: PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,optionname,ifields,&nfields,&flg);
294: PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,optionname_col,ifields_col,&nfields_col,&flg_col);
295: if (!flg) break;
296: else if (flg && !flg_col) {
297: if (!nfields) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields");
298: PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields);
299: } else {
300: if (!nfields || !nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields");
301: if (nfields != nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Number of row and column fields must match");
302: PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields_col);
303: }
304: }
305: if (i > 0) {
306: /* Makes command-line setting of splits take precedence over setting them in code.
307: Otherwise subsequent calls to PCFieldSplitSetIS() or PCFieldSplitSetFields() would
308: create new splits, which would probably not be what the user wanted. */
309: jac->splitdefined = PETSC_TRUE;
310: }
311: PetscFree(ifields);
312: PetscFree(ifields_col);
313: return(0);
314: }
316: static PetscErrorCode PCFieldSplitSetDefaults(PC pc)
317: {
318: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
319: PetscErrorCode ierr;
320: PC_FieldSplitLink ilink = jac->head;
321: PetscBool fieldsplit_default = PETSC_FALSE,stokes = PETSC_FALSE,coupling = PETSC_FALSE;
322: PetscInt i;
325: /*
326: Kinda messy, but at least this now uses DMCreateFieldDecomposition().
327: Should probably be rewritten.
328: */
329: if (!ilink) {
330: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);
331: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_coupling",&coupling,NULL);
332: if (pc->dm && jac->dm_splits && !stokes && !coupling) {
333: PetscInt numFields, f, i, j;
334: char **fieldNames;
335: IS *fields;
336: DM *dms;
337: DM subdm[128];
338: PetscBool flg;
340: DMCreateFieldDecomposition(pc->dm, &numFields, &fieldNames, &fields, &dms);
341: /* Allow the user to prescribe the splits */
342: for (i = 0, flg = PETSC_TRUE;; i++) {
343: PetscInt ifields[128];
344: IS compField;
345: char optionname[128], splitname[8];
346: PetscInt nfields = numFields;
348: PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%D_fields", i);
349: PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix, optionname, ifields, &nfields, &flg);
350: if (!flg) break;
351: if (numFields > 128) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot currently support %d > 128 fields", numFields);
352: DMCreateSubDM(pc->dm, nfields, ifields, &compField, &subdm[i]);
353: if (nfields == 1) {
354: PCFieldSplitSetIS(pc, fieldNames[ifields[0]], compField);
355: } else {
356: PetscSNPrintf(splitname, sizeof(splitname), "%D", i);
357: PCFieldSplitSetIS(pc, splitname, compField);
358: }
359: ISDestroy(&compField);
360: for (j = 0; j < nfields; ++j) {
361: f = ifields[j];
362: PetscFree(fieldNames[f]);
363: ISDestroy(&fields[f]);
364: }
365: }
366: if (i == 0) {
367: for (f = 0; f < numFields; ++f) {
368: PCFieldSplitSetIS(pc, fieldNames[f], fields[f]);
369: PetscFree(fieldNames[f]);
370: ISDestroy(&fields[f]);
371: }
372: } else {
373: for (j=0; j<numFields; j++) {
374: DMDestroy(dms+j);
375: }
376: PetscFree(dms);
377: PetscMalloc1(i, &dms);
378: for (j = 0; j < i; ++j) dms[j] = subdm[j];
379: }
380: PetscFree(fieldNames);
381: PetscFree(fields);
382: if (dms) {
383: PetscInfo(pc, "Setting up physics based fieldsplit preconditioner using the embedded DM\n");
384: for (ilink = jac->head, i = 0; ilink; ilink = ilink->next, ++i) {
385: const char *prefix;
386: PetscObjectGetOptionsPrefix((PetscObject)(ilink->ksp),&prefix);
387: PetscObjectSetOptionsPrefix((PetscObject)(dms[i]), prefix);
388: KSPSetDM(ilink->ksp, dms[i]);
389: KSPSetDMActive(ilink->ksp, PETSC_FALSE);
390: PetscObjectIncrementTabLevel((PetscObject)dms[i],(PetscObject)ilink->ksp,0);
391: DMDestroy(&dms[i]);
392: }
393: PetscFree(dms);
394: }
395: } else {
396: if (jac->bs <= 0) {
397: if (pc->pmat) {
398: MatGetBlockSize(pc->pmat,&jac->bs);
399: } else jac->bs = 1;
400: }
402: if (stokes) {
403: IS zerodiags,rest;
404: PetscInt nmin,nmax;
406: MatGetOwnershipRange(pc->mat,&nmin,&nmax);
407: MatFindZeroDiagonals(pc->mat,&zerodiags);
408: ISComplement(zerodiags,nmin,nmax,&rest);
409: PCFieldSplitSetIS(pc,"0",rest);
410: PCFieldSplitSetIS(pc,"1",zerodiags);
411: ISDestroy(&zerodiags);
412: ISDestroy(&rest);
413: } else if (coupling) {
414: IS coupling,rest;
415: PetscInt nmin,nmax;
417: MatGetOwnershipRange(pc->mat,&nmin,&nmax);
418: MatFindOffBlockDiagonalEntries(pc->mat,&coupling);
419: ISCreateStride(PetscObjectComm((PetscObject)pc->mat),nmax-nmin,nmin,1,&rest);
420: ISSetIdentity(rest);
421: PCFieldSplitSetIS(pc,"0",rest);
422: PCFieldSplitSetIS(pc,"1",coupling);
423: ISDestroy(&coupling);
424: ISDestroy(&rest);
425: } else {
426: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&fieldsplit_default,NULL);
427: if (!fieldsplit_default) {
428: /* Allow user to set fields from command line, if bs was known at the time of PCSetFromOptions_FieldSplit()
429: then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */
430: PCFieldSplitSetRuntimeSplits_Private(pc);
431: if (jac->splitdefined) {PetscInfo(pc,"Splits defined using the options database\n");}
432: }
433: if ((fieldsplit_default || !jac->splitdefined) && !jac->isrestrict) {
434: PetscInfo(pc,"Using default splitting of fields\n");
435: for (i=0; i<jac->bs; i++) {
436: char splitname[8];
437: PetscSNPrintf(splitname,sizeof(splitname),"%D",i);
438: PCFieldSplitSetFields(pc,splitname,1,&i,&i);
439: }
440: jac->defaultsplit = PETSC_TRUE;
441: }
442: }
443: }
444: } else if (jac->nsplits == 1) {
445: if (ilink->is) {
446: IS is2;
447: PetscInt nmin,nmax;
449: MatGetOwnershipRange(pc->mat,&nmin,&nmax);
450: ISComplement(ilink->is,nmin,nmax,&is2);
451: PCFieldSplitSetIS(pc,"1",is2);
452: ISDestroy(&is2);
453: } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Must provide at least two sets of fields to PCFieldSplit()");
454: }
456: if (jac->nsplits < 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Unhandled case, must have at least two fields, not %d", jac->nsplits);
457: return(0);
458: }
460: PETSC_EXTERN PetscErrorCode PetscOptionsFindPairPrefix_Private(PetscOptions,const char pre[], const char name[], char *value[], PetscBool *flg);
462: static PetscErrorCode PCSetUp_FieldSplit(PC pc)
463: {
464: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
465: PetscErrorCode ierr;
466: PC_FieldSplitLink ilink;
467: PetscInt i,nsplit;
468: PetscBool sorted, sorted_col;
471: PCFieldSplitSetDefaults(pc);
472: nsplit = jac->nsplits;
473: ilink = jac->head;
475: /* get the matrices for each split */
476: if (!jac->issetup) {
477: PetscInt rstart,rend,nslots,bs;
479: jac->issetup = PETSC_TRUE;
481: /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */
482: if (jac->defaultsplit || !ilink->is) {
483: if (jac->bs <= 0) jac->bs = nsplit;
484: }
485: bs = jac->bs;
486: MatGetOwnershipRange(pc->pmat,&rstart,&rend);
487: nslots = (rend - rstart)/bs;
488: for (i=0; i<nsplit; i++) {
489: if (jac->defaultsplit) {
490: ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+i,nsplit,&ilink->is);
491: ISDuplicate(ilink->is,&ilink->is_col);
492: } else if (!ilink->is) {
493: if (ilink->nfields > 1) {
494: PetscInt *ii,*jj,j,k,nfields = ilink->nfields,*fields = ilink->fields,*fields_col = ilink->fields_col;
495: PetscMalloc1(ilink->nfields*nslots,&ii);
496: PetscMalloc1(ilink->nfields*nslots,&jj);
497: for (j=0; j<nslots; j++) {
498: for (k=0; k<nfields; k++) {
499: ii[nfields*j + k] = rstart + bs*j + fields[k];
500: jj[nfields*j + k] = rstart + bs*j + fields_col[k];
501: }
502: }
503: ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,ii,PETSC_OWN_POINTER,&ilink->is);
504: ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,jj,PETSC_OWN_POINTER,&ilink->is_col);
505: ISSetBlockSize(ilink->is, nfields);
506: ISSetBlockSize(ilink->is_col, nfields);
507: } else {
508: ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields[0],bs,&ilink->is);
509: ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields_col[0],bs,&ilink->is_col);
510: }
511: }
512: ISSorted(ilink->is,&sorted);
513: if (ilink->is_col) { ISSorted(ilink->is_col,&sorted_col); }
514: if (!sorted || !sorted_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Fields must be sorted when creating split");
515: ilink = ilink->next;
516: }
517: }
519: ilink = jac->head;
520: if (!jac->pmat) {
521: Vec xtmp;
523: MatCreateVecs(pc->pmat,&xtmp,NULL);
524: PetscMalloc1(nsplit,&jac->pmat);
525: PetscMalloc2(nsplit,&jac->x,nsplit,&jac->y);
526: for (i=0; i<nsplit; i++) {
527: MatNullSpace sp;
529: /* Check for preconditioning matrix attached to IS */
530: PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &jac->pmat[i]);
531: if (jac->pmat[i]) {
532: PetscObjectReference((PetscObject) jac->pmat[i]);
533: if (jac->type == PC_COMPOSITE_SCHUR) {
534: jac->schur_user = jac->pmat[i];
536: PetscObjectReference((PetscObject) jac->schur_user);
537: }
538: } else {
539: const char *prefix;
540: MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->pmat[i]);
541: KSPGetOptionsPrefix(ilink->ksp,&prefix);
542: MatSetOptionsPrefix(jac->pmat[i],prefix);
543: MatViewFromOptions(jac->pmat[i],NULL,"-mat_view");
544: }
545: /* create work vectors for each split */
546: MatCreateVecs(jac->pmat[i],&jac->x[i],&jac->y[i]);
547: ilink->x = jac->x[i]; ilink->y = jac->y[i]; ilink->z = NULL;
548: /* compute scatter contexts needed by multiplicative versions and non-default splits */
549: VecScatterCreate(xtmp,ilink->is,jac->x[i],NULL,&ilink->sctx);
550: PetscObjectQuery((PetscObject) ilink->is, "nearnullspace", (PetscObject*) &sp);
551: if (sp) {
552: MatSetNearNullSpace(jac->pmat[i], sp);
553: }
554: ilink = ilink->next;
555: }
556: VecDestroy(&xtmp);
557: } else {
558: for (i=0; i<nsplit; i++) {
559: Mat pmat;
561: /* Check for preconditioning matrix attached to IS */
562: PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &pmat);
563: if (!pmat) {
564: MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_REUSE_MATRIX,&jac->pmat[i]);
565: }
566: ilink = ilink->next;
567: }
568: }
569: if (jac->diag_use_amat) {
570: ilink = jac->head;
571: if (!jac->mat) {
572: PetscMalloc1(nsplit,&jac->mat);
573: for (i=0; i<nsplit; i++) {
574: MatCreateSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->mat[i]);
575: ilink = ilink->next;
576: }
577: } else {
578: for (i=0; i<nsplit; i++) {
579: if (jac->mat[i]) {MatCreateSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_REUSE_MATRIX,&jac->mat[i]);}
580: ilink = ilink->next;
581: }
582: }
583: } else {
584: jac->mat = jac->pmat;
585: }
587: /* Check for null space attached to IS */
588: ilink = jac->head;
589: for (i=0; i<nsplit; i++) {
590: MatNullSpace sp;
592: PetscObjectQuery((PetscObject) ilink->is, "nullspace", (PetscObject*) &sp);
593: if (sp) {
594: MatSetNullSpace(jac->mat[i], sp);
595: }
596: ilink = ilink->next;
597: }
599: if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR) {
600: /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */
601: /* FIXME: Can/should we reuse jac->mat whenever (jac->diag_use_amat) is true? */
602: ilink = jac->head;
603: if (nsplit == 2 && jac->type == PC_COMPOSITE_MULTIPLICATIVE) {
604: /* special case need where Afield[0] is not needed and only certain columns of Afield[1] are needed since update is only on those rows of the solution */
605: if (!jac->Afield) {
606: PetscCalloc1(nsplit,&jac->Afield);
607: if (jac->offdiag_use_amat) {
608: MatCreateSubMatrix(pc->mat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);
609: } else {
610: MatCreateSubMatrix(pc->pmat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);
611: }
612: } else {
613: if (jac->offdiag_use_amat) {
614: MatCreateSubMatrix(pc->mat,ilink->next->is,ilink->is,MAT_REUSE_MATRIX,&jac->Afield[1]);
615: } else {
616: MatCreateSubMatrix(pc->pmat,ilink->next->is,ilink->is,MAT_REUSE_MATRIX,&jac->Afield[1]);
617: }
618: }
619: } else {
620: if (!jac->Afield) {
621: PetscMalloc1(nsplit,&jac->Afield);
622: for (i=0; i<nsplit; i++) {
623: if (jac->offdiag_use_amat) {
624: MatCreateSubMatrix(pc->mat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);
625: } else {
626: MatCreateSubMatrix(pc->pmat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);
627: }
628: ilink = ilink->next;
629: }
630: } else {
631: for (i=0; i<nsplit; i++) {
632: if (jac->offdiag_use_amat) {
633: MatCreateSubMatrix(pc->mat,ilink->is,NULL,MAT_REUSE_MATRIX,&jac->Afield[i]);
634: } else {
635: MatCreateSubMatrix(pc->pmat,ilink->is,NULL,MAT_REUSE_MATRIX,&jac->Afield[i]);
636: }
637: ilink = ilink->next;
638: }
639: }
640: }
641: }
643: if (jac->type == PC_COMPOSITE_SCHUR) {
644: IS ccis;
645: PetscBool isspd;
646: PetscInt rstart,rend;
647: char lscname[256];
648: PetscObject LSC_L;
650: if (nsplit != 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"To use Schur complement preconditioner you must have exactly 2 fields");
652: /* If pc->mat is SPD, don't scale by -1 the Schur complement */
653: if (jac->schurscale == (PetscScalar)-1.0) {
654: MatGetOption(pc->pmat,MAT_SPD,&isspd);
655: jac->schurscale = (isspd == PETSC_TRUE) ? 1.0 : -1.0;
656: }
658: /* When extracting off-diagonal submatrices, we take complements from this range */
659: MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);
661: /* need to handle case when one is resetting up the preconditioner */
662: if (jac->schur) {
663: KSP kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper;
665: MatSchurComplementGetKSP(jac->schur, &kspInner);
666: ilink = jac->head;
667: ISComplement(ilink->is_col,rstart,rend,&ccis);
668: if (jac->offdiag_use_amat) {
669: MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);
670: } else {
671: MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);
672: }
673: ISDestroy(&ccis);
674: ilink = ilink->next;
675: ISComplement(ilink->is_col,rstart,rend,&ccis);
676: if (jac->offdiag_use_amat) {
677: MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);
678: } else {
679: MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);
680: }
681: ISDestroy(&ccis);
682: MatSchurComplementUpdateSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);
683: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
684: MatDestroy(&jac->schurp);
685: MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);
686: }
687: if (kspA != kspInner) {
688: KSPSetOperators(kspA,jac->mat[0],jac->pmat[0]);
689: }
690: if (kspUpper != kspA) {
691: KSPSetOperators(kspUpper,jac->mat[0],jac->pmat[0]);
692: }
693: KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));
694: } else {
695: const char *Dprefix;
696: char schurprefix[256], schurmatprefix[256];
697: char schurtestoption[256];
698: MatNullSpace sp;
699: PetscBool flg;
701: /* extract the A01 and A10 matrices */
702: ilink = jac->head;
703: ISComplement(ilink->is_col,rstart,rend,&ccis);
704: if (jac->offdiag_use_amat) {
705: MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
706: } else {
707: MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
708: }
709: ISDestroy(&ccis);
710: ilink = ilink->next;
711: ISComplement(ilink->is_col,rstart,rend,&ccis);
712: if (jac->offdiag_use_amat) {
713: MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
714: } else {
715: MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
716: }
717: ISDestroy(&ccis);
719: /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */
720: MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur);
721: MatSetType(jac->schur,MATSCHURCOMPLEMENT);
722: MatSchurComplementSetSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);
723: PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
724: /* Note that the inner KSP is NOT going to inherit this prefix, and if it did, it would be reset just below. Is that what we want? */
725: MatSetOptionsPrefix(jac->schur,schurmatprefix);
726: MatSetFromOptions(jac->schur);
727: MatGetNullSpace(jac->mat[1], &sp);
728: if (sp) {
729: MatSetNullSpace(jac->schur, sp);
730: }
732: PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);
733: PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
734: if (flg) {
735: DM dmInner;
736: KSP kspInner;
738: MatSchurComplementGetKSP(jac->schur, &kspInner);
739: PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
740: /* Indent this deeper to emphasize the "inner" nature of this solver. */
741: PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject) pc, 2);
742: KSPSetOptionsPrefix(kspInner, schurprefix);
744: /* Set DM for new solver */
745: KSPGetDM(jac->head->ksp, &dmInner);
746: KSPSetDM(kspInner, dmInner);
747: KSPSetDMActive(kspInner, PETSC_FALSE);
748: } else {
749: /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or
750: * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact,
751: * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for
752: * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make
753: * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used
754: * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */
755: KSPSetType(jac->head->ksp,KSPGMRES);
756: MatSchurComplementSetKSP(jac->schur,jac->head->ksp);
757: }
758: KSPSetOperators(jac->head->ksp,jac->mat[0],jac->pmat[0]);
759: KSPSetFromOptions(jac->head->ksp);
760: MatSetFromOptions(jac->schur);
762: PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);
763: PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
764: if (flg) {
765: DM dmInner;
767: PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
768: KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper);
769: KSPSetErrorIfNotConverged(jac->kspupper,pc->erroriffailure);
770: KSPSetOptionsPrefix(jac->kspupper, schurprefix);
771: KSPGetDM(jac->head->ksp, &dmInner);
772: KSPSetDM(jac->kspupper, dmInner);
773: KSPSetDMActive(jac->kspupper, PETSC_FALSE);
774: KSPSetFromOptions(jac->kspupper);
775: KSPSetOperators(jac->kspupper,jac->mat[0],jac->pmat[0]);
776: VecDuplicate(jac->head->x, &jac->head->z);
777: } else {
778: jac->kspupper = jac->head->ksp;
779: PetscObjectReference((PetscObject) jac->head->ksp);
780: }
782: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
783: MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);
784: }
785: KSPCreate(PetscObjectComm((PetscObject)pc),&jac->kspschur);
786: KSPSetErrorIfNotConverged(jac->kspschur,pc->erroriffailure);
787: PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);
788: PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);
789: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) {
790: PC pcschur;
791: KSPGetPC(jac->kspschur,&pcschur);
792: PCSetType(pcschur,PCNONE);
793: /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */
794: } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) {
795: MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user);
796: }
797: KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));
798: KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix);
799: KSPSetOptionsPrefix(jac->kspschur, Dprefix);
800: /* propagate DM */
801: {
802: DM sdm;
803: KSPGetDM(jac->head->next->ksp, &sdm);
804: if (sdm) {
805: KSPSetDM(jac->kspschur, sdm);
806: KSPSetDMActive(jac->kspschur, PETSC_FALSE);
807: }
808: }
809: /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
810: /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */
811: KSPSetFromOptions(jac->kspschur);
812: }
814: /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */
815: PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_L",ilink->splitname);
816: PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);
817: if (!LSC_L) {PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);}
818: if (LSC_L) {PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L);}
819: PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_Lp",ilink->splitname);
820: PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);
821: if (!LSC_L) {PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);}
822: if (LSC_L) {PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L);}
823: } else {
824: /* set up the individual splits' PCs */
825: i = 0;
826: ilink = jac->head;
827: while (ilink) {
828: KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i]);
829: /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
830: if (!jac->suboptionsset) {KSPSetFromOptions(ilink->ksp);}
831: i++;
832: ilink = ilink->next;
833: }
834: }
836: jac->suboptionsset = PETSC_TRUE;
837: return(0);
838: }
840: #define FieldSplitSplitSolveAdd(ilink,xx,yy) \
841: (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
842: VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
843: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\
844: KSPSolve(ilink->ksp,ilink->x,ilink->y) || \
845: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\
846: VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \
847: VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE))
849: static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y)
850: {
851: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
852: PetscErrorCode ierr;
853: PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
854: KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;
857: switch (jac->schurfactorization) {
858: case PC_FIELDSPLIT_SCHUR_FACT_DIAG:
859: /* [A00 0; 0 -S], positive definite, suitable for MINRES */
860: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
861: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
862: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
863: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
864: KSPSolve(kspA,ilinkA->x,ilinkA->y);
865: PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
866: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
867: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
868: PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
869: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
870: PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
871: VecScale(ilinkD->y,jac->schurscale);
872: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
873: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
874: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
875: break;
876: case PC_FIELDSPLIT_SCHUR_FACT_LOWER:
877: /* [A00 0; A10 S], suitable for left preconditioning */
878: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
879: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
880: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
881: KSPSolve(kspA,ilinkA->x,ilinkA->y);
882: PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
883: MatMult(jac->C,ilinkA->y,ilinkD->x);
884: VecScale(ilinkD->x,-1.);
885: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
886: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
887: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
888: PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
889: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
890: PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
891: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
892: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
893: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
894: break;
895: case PC_FIELDSPLIT_SCHUR_FACT_UPPER:
896: /* [A00 A01; 0 S], suitable for right preconditioning */
897: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
898: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
899: PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
900: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
901: PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL); MatMult(jac->B,ilinkD->y,ilinkA->x);
902: VecScale(ilinkA->x,-1.);
903: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);
904: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
905: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);
906: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
907: KSPSolve(kspA,ilinkA->x,ilinkA->y);
908: PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
909: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
910: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
911: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
912: break;
913: case PC_FIELDSPLIT_SCHUR_FACT_FULL:
914: /* [1 0; A10 A00^{-1} 1] [A00 0; 0 S] [1 A00^{-1}A01; 0 1], an exact solve if applied exactly, needs one extra solve with A */
915: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
916: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
917: PetscLogEventBegin(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);
918: KSPSolve(kspLower,ilinkA->x,ilinkA->y);
919: PetscLogEventEnd(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);
920: MatMult(jac->C,ilinkA->y,ilinkD->x);
921: VecScale(ilinkD->x,-1.0);
922: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
923: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
925: PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
926: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
927: PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
928: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
929: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
931: if (kspUpper == kspA) {
932: MatMult(jac->B,ilinkD->y,ilinkA->y);
933: VecAXPY(ilinkA->x,-1.0,ilinkA->y);
934: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
935: KSPSolve(kspA,ilinkA->x,ilinkA->y);
936: PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
937: } else {
938: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
939: KSPSolve(kspA,ilinkA->x,ilinkA->y);
940: PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
941: MatMult(jac->B,ilinkD->y,ilinkA->x);
942: PetscLogEventBegin(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);
943: KSPSolve(kspUpper,ilinkA->x,ilinkA->z);
944: PetscLogEventEnd(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);
945: VecAXPY(ilinkA->y,-1.0,ilinkA->z);
946: }
947: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
948: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
949: }
950: return(0);
951: }
953: static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y)
954: {
955: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
956: PetscErrorCode ierr;
957: PC_FieldSplitLink ilink = jac->head;
958: PetscInt cnt,bs;
959: KSPConvergedReason reason;
962: if (jac->type == PC_COMPOSITE_ADDITIVE) {
963: if (jac->defaultsplit) {
964: VecGetBlockSize(x,&bs);
965: if (jac->bs > 0 && bs != jac->bs) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Blocksize of x vector %D does not match fieldsplit blocksize %D",bs,jac->bs);
966: VecGetBlockSize(y,&bs);
967: if (jac->bs > 0 && bs != jac->bs) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Blocksize of y vector %D does not match fieldsplit blocksize %D",bs,jac->bs);
968: VecStrideGatherAll(x,jac->x,INSERT_VALUES);
969: while (ilink) {
970: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
971: KSPSolve(ilink->ksp,ilink->x,ilink->y);
972: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
973: KSPGetConvergedReason(ilink->ksp,&reason);
974: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
975: pc->failedreason = PC_SUBPC_ERROR;
976: }
977: ilink = ilink->next;
978: }
979: VecStrideScatterAll(jac->y,y,INSERT_VALUES);
980: } else {
981: VecSet(y,0.0);
982: while (ilink) {
983: FieldSplitSplitSolveAdd(ilink,x,y);
984: KSPGetConvergedReason(ilink->ksp,&reason);
985: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
986: pc->failedreason = PC_SUBPC_ERROR;
987: }
988: ilink = ilink->next;
989: }
990: }
991: } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) {
992: VecSet(y,0.0);
993: /* solve on first block for first block variables */
994: VecScatterBegin(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);
995: VecScatterEnd(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);
996: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
997: KSPSolve(ilink->ksp,ilink->x,ilink->y);
998: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
999: KSPGetConvergedReason(ilink->ksp,&reason);
1000: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1001: pc->failedreason = PC_SUBPC_ERROR;
1002: }
1003: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1004: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1006: /* compute the residual only onto second block variables using first block variables */
1007: MatMult(jac->Afield[1],ilink->y,ilink->next->x);
1008: ilink = ilink->next;
1009: VecScale(ilink->x,-1.0);
1010: VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1011: VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1013: /* solve on second block variables */
1014: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1015: KSPSolve(ilink->ksp,ilink->x,ilink->y);
1016: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1017: KSPGetConvergedReason(ilink->ksp,&reason);
1018: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1019: pc->failedreason = PC_SUBPC_ERROR;
1020: }
1021: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1022: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1023: } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1024: if (!jac->w1) {
1025: VecDuplicate(x,&jac->w1);
1026: VecDuplicate(x,&jac->w2);
1027: }
1028: VecSet(y,0.0);
1029: FieldSplitSplitSolveAdd(ilink,x,y);
1030: KSPGetConvergedReason(ilink->ksp,&reason);
1031: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1032: pc->failedreason = PC_SUBPC_ERROR;
1033: }
1034: cnt = 1;
1035: while (ilink->next) {
1036: ilink = ilink->next;
1037: /* compute the residual only over the part of the vector needed */
1038: MatMult(jac->Afield[cnt++],y,ilink->x);
1039: VecScale(ilink->x,-1.0);
1040: VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1041: VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1042: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1043: KSPSolve(ilink->ksp,ilink->x,ilink->y);
1044: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1045: KSPGetConvergedReason(ilink->ksp,&reason);
1046: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1047: pc->failedreason = PC_SUBPC_ERROR;
1048: }
1049: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1050: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1051: }
1052: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1053: cnt -= 2;
1054: while (ilink->previous) {
1055: ilink = ilink->previous;
1056: /* compute the residual only over the part of the vector needed */
1057: MatMult(jac->Afield[cnt--],y,ilink->x);
1058: VecScale(ilink->x,-1.0);
1059: VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1060: VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1061: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1062: KSPSolve(ilink->ksp,ilink->x,ilink->y);
1063: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1064: KSPGetConvergedReason(ilink->ksp,&reason);
1065: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1066: pc->failedreason = PC_SUBPC_ERROR;
1067: }
1068: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1069: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1070: }
1071: }
1072: } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type);
1073: return(0);
1074: }
1076: #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \
1077: (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
1078: VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
1079: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \
1080: KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \
1081: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \
1082: VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \
1083: VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE))
1085: static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y)
1086: {
1087: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1088: PetscErrorCode ierr;
1089: PC_FieldSplitLink ilink = jac->head;
1090: PetscInt bs;
1091: KSPConvergedReason reason;
1094: if (jac->type == PC_COMPOSITE_ADDITIVE) {
1095: if (jac->defaultsplit) {
1096: VecGetBlockSize(x,&bs);
1097: if (jac->bs > 0 && bs != jac->bs) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Blocksize of x vector %D does not match fieldsplit blocksize %D",bs,jac->bs);
1098: VecGetBlockSize(y,&bs);
1099: if (jac->bs > 0 && bs != jac->bs) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Blocksize of y vector %D does not match fieldsplit blocksize %D",bs,jac->bs);
1100: VecStrideGatherAll(x,jac->x,INSERT_VALUES);
1101: while (ilink) {
1102: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1103: KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);
1104: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1105: KSPGetConvergedReason(ilink->ksp,&reason);
1106: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1107: pc->failedreason = PC_SUBPC_ERROR;
1108: }
1109: ilink = ilink->next;
1110: }
1111: VecStrideScatterAll(jac->y,y,INSERT_VALUES);
1112: } else {
1113: VecSet(y,0.0);
1114: while (ilink) {
1115: FieldSplitSplitSolveAddTranspose(ilink,x,y);
1116: KSPGetConvergedReason(ilink->ksp,&reason);
1117: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1118: pc->failedreason = PC_SUBPC_ERROR;
1119: }
1120: ilink = ilink->next;
1121: }
1122: }
1123: } else {
1124: if (!jac->w1) {
1125: VecDuplicate(x,&jac->w1);
1126: VecDuplicate(x,&jac->w2);
1127: }
1128: VecSet(y,0.0);
1129: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1130: FieldSplitSplitSolveAddTranspose(ilink,x,y);
1131: KSPGetConvergedReason(ilink->ksp,&reason);
1132: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1133: pc->failedreason = PC_SUBPC_ERROR;
1134: }
1135: while (ilink->next) {
1136: ilink = ilink->next;
1137: MatMultTranspose(pc->mat,y,jac->w1);
1138: VecWAXPY(jac->w2,-1.0,jac->w1,x);
1139: FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1140: }
1141: while (ilink->previous) {
1142: ilink = ilink->previous;
1143: MatMultTranspose(pc->mat,y,jac->w1);
1144: VecWAXPY(jac->w2,-1.0,jac->w1,x);
1145: FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1146: }
1147: } else {
1148: while (ilink->next) { /* get to last entry in linked list */
1149: ilink = ilink->next;
1150: }
1151: FieldSplitSplitSolveAddTranspose(ilink,x,y);
1152: KSPGetConvergedReason(ilink->ksp,&reason);
1153: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1154: pc->failedreason = PC_SUBPC_ERROR;
1155: }
1156: while (ilink->previous) {
1157: ilink = ilink->previous;
1158: MatMultTranspose(pc->mat,y,jac->w1);
1159: VecWAXPY(jac->w2,-1.0,jac->w1,x);
1160: FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1161: }
1162: }
1163: }
1164: return(0);
1165: }
1167: static PetscErrorCode PCReset_FieldSplit(PC pc)
1168: {
1169: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1170: PetscErrorCode ierr;
1171: PC_FieldSplitLink ilink = jac->head,next;
1174: while (ilink) {
1175: KSPDestroy(&ilink->ksp);
1176: VecDestroy(&ilink->x);
1177: VecDestroy(&ilink->y);
1178: VecDestroy(&ilink->z);
1179: VecScatterDestroy(&ilink->sctx);
1180: ISDestroy(&ilink->is);
1181: ISDestroy(&ilink->is_col);
1182: PetscFree(ilink->splitname);
1183: PetscFree(ilink->fields);
1184: PetscFree(ilink->fields_col);
1185: next = ilink->next;
1186: PetscFree(ilink);
1187: ilink = next;
1188: }
1189: jac->head = NULL;
1190: PetscFree2(jac->x,jac->y);
1191: if (jac->mat && jac->mat != jac->pmat) {
1192: MatDestroyMatrices(jac->nsplits,&jac->mat);
1193: } else if (jac->mat) {
1194: jac->mat = NULL;
1195: }
1196: if (jac->pmat) {MatDestroyMatrices(jac->nsplits,&jac->pmat);}
1197: if (jac->Afield) {MatDestroyMatrices(jac->nsplits,&jac->Afield);}
1198: jac->nsplits = 0;
1199: VecDestroy(&jac->w1);
1200: VecDestroy(&jac->w2);
1201: MatDestroy(&jac->schur);
1202: MatDestroy(&jac->schurp);
1203: MatDestroy(&jac->schur_user);
1204: KSPDestroy(&jac->kspschur);
1205: KSPDestroy(&jac->kspupper);
1206: MatDestroy(&jac->B);
1207: MatDestroy(&jac->C);
1208: jac->isrestrict = PETSC_FALSE;
1209: return(0);
1210: }
1212: static PetscErrorCode PCDestroy_FieldSplit(PC pc)
1213: {
1214: PetscErrorCode ierr;
1217: PCReset_FieldSplit(pc);
1218: PetscFree(pc->data);
1219: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL);
1220: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL);
1221: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL);
1222: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL);
1223: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL);
1224: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL);
1225: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL);
1226: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL);
1227: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",NULL);
1228: return(0);
1229: }
1231: static PetscErrorCode PCSetFromOptions_FieldSplit(PetscOptionItems *PetscOptionsObject,PC pc)
1232: {
1233: PetscErrorCode ierr;
1234: PetscInt bs;
1235: PetscBool flg,stokes = PETSC_FALSE;
1236: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1237: PCCompositeType ctype;
1240: PetscOptionsHead(PetscOptionsObject,"FieldSplit options");
1241: PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL);
1242: PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);
1243: if (flg) {
1244: PCFieldSplitSetBlockSize(pc,bs);
1245: }
1246: jac->diag_use_amat = pc->useAmat;
1247: PetscOptionsBool("-pc_fieldsplit_diag_use_amat","Use Amat (not Pmat) to extract diagonal fieldsplit blocks", "PCFieldSplitSetDiagUseAmat",jac->diag_use_amat,&jac->diag_use_amat,NULL);
1248: jac->offdiag_use_amat = pc->useAmat;
1249: PetscOptionsBool("-pc_fieldsplit_off_diag_use_amat","Use Amat (not Pmat) to extract off-diagonal fieldsplit blocks", "PCFieldSplitSetOffDiagUseAmat",jac->offdiag_use_amat,&jac->offdiag_use_amat,NULL);
1250: /* FIXME: No programmatic equivalent to the following. */
1251: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);
1252: if (stokes) {
1253: PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);
1254: jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF;
1255: }
1257: PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);
1258: if (flg) {
1259: PCFieldSplitSetType(pc,ctype);
1260: }
1261: /* Only setup fields once */
1262: if ((jac->bs > 0) && (jac->nsplits == 0)) {
1263: /* only allow user to set fields from command line if bs is already known.
1264: otherwise user can set them in PCFieldSplitSetDefaults() */
1265: PCFieldSplitSetRuntimeSplits_Private(pc);
1266: if (jac->splitdefined) {PetscInfo(pc,"Splits defined using the options database\n");}
1267: }
1268: if (jac->type == PC_COMPOSITE_SCHUR) {
1269: PetscOptionsGetEnum(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);
1270: if (flg) {PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");}
1271: PetscOptionsEnum("-pc_fieldsplit_schur_fact_type","Which off-diagonal parts of the block factorization to use","PCFieldSplitSetSchurFactType",PCFieldSplitSchurFactTypes,(PetscEnum)jac->schurfactorization,(PetscEnum*)&jac->schurfactorization,NULL);
1272: PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSetSchurPre",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL);
1273: PetscOptionsScalar("-pc_fieldsplit_schur_scale","Scale Schur complement","PCFieldSplitSetSchurScale",jac->schurscale,&jac->schurscale,NULL);
1274: }
1275: PetscOptionsTail();
1276: return(0);
1277: }
1279: /*------------------------------------------------------------------------------------*/
1281: static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1282: {
1283: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1284: PetscErrorCode ierr;
1285: PC_FieldSplitLink ilink,next = jac->head;
1286: char prefix[128];
1287: PetscInt i;
1290: if (jac->splitdefined) {
1291: PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1292: return(0);
1293: }
1294: for (i=0; i<n; i++) {
1295: if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs);
1296: if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]);
1297: }
1298: PetscNew(&ilink);
1299: if (splitname) {
1300: PetscStrallocpy(splitname,&ilink->splitname);
1301: } else {
1302: PetscMalloc1(3,&ilink->splitname);
1303: PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);
1304: }
1305: ilink->event = jac->nsplits < 5 ? KSP_Solve_FS_0 + jac->nsplits : KSP_Solve_FS_0 + 4; /* Any split great than 4 gets logged in the 4th split */
1306: PetscMalloc1(n,&ilink->fields);
1307: PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));
1308: PetscMalloc1(n,&ilink->fields_col);
1309: PetscMemcpy(ilink->fields_col,fields_col,n*sizeof(PetscInt));
1311: ilink->nfields = n;
1312: ilink->next = NULL;
1313: KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1314: KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);
1315: PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1316: KSPSetType(ilink->ksp,KSPPREONLY);
1317: PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);
1319: PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);
1320: KSPSetOptionsPrefix(ilink->ksp,prefix);
1322: if (!next) {
1323: jac->head = ilink;
1324: ilink->previous = NULL;
1325: } else {
1326: while (next->next) {
1327: next = next->next;
1328: }
1329: next->next = ilink;
1330: ilink->previous = next;
1331: }
1332: jac->nsplits++;
1333: return(0);
1334: }
1336: static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp)
1337: {
1338: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1342: if (!jac->schur) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitGetSubKSP()");
1343: PetscMalloc1(jac->nsplits,subksp);
1344: MatSchurComplementGetKSP(jac->schur,*subksp);
1346: (*subksp)[1] = jac->kspschur;
1347: if (n) *n = jac->nsplits;
1348: return(0);
1349: }
1351: static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp)
1352: {
1353: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1354: PetscErrorCode ierr;
1355: PetscInt cnt = 0;
1356: PC_FieldSplitLink ilink = jac->head;
1359: PetscMalloc1(jac->nsplits,subksp);
1360: while (ilink) {
1361: (*subksp)[cnt++] = ilink->ksp;
1362: ilink = ilink->next;
1363: }
1364: if (cnt != jac->nsplits) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupt PCFIELDSPLIT object: number of splits in linked list %D does not match number in object %D",cnt,jac->nsplits);
1365: if (n) *n = jac->nsplits;
1366: return(0);
1367: }
1369: /*@C
1370: PCFieldSplitRestrictIS - Restricts the fieldsplit ISs to be within a given IS.
1372: Input Parameters:
1373: + pc - the preconditioner context
1374: + is - the index set that defines the indices to which the fieldsplit is to be restricted
1376: Level: advanced
1378: @*/
1379: PetscErrorCode PCFieldSplitRestrictIS(PC pc,IS isy)
1380: {
1386: PetscTryMethod(pc,"PCFieldSplitRestrictIS_C",(PC,IS),(pc,isy));
1387: return(0);
1388: }
1391: static PetscErrorCode PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy)
1392: {
1393: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1394: PetscErrorCode ierr;
1395: PC_FieldSplitLink ilink = jac->head, next;
1396: PetscInt localsize,size,sizez,i;
1397: const PetscInt *ind, *indz;
1398: PetscInt *indc, *indcz;
1399: PetscBool flg;
1402: ISGetLocalSize(isy,&localsize);
1403: MPI_Scan(&localsize,&size,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)isy));
1404: size -= localsize;
1405: while(ilink) {
1406: IS isrl,isr;
1407: PC subpc;
1408: ISEmbed(ilink->is, isy, PETSC_TRUE, &isrl);
1409: ISGetLocalSize(isrl,&localsize);
1410: PetscMalloc1(localsize,&indc);
1411: ISGetIndices(isrl,&ind);
1412: PetscMemcpy(indc,ind,localsize*sizeof(PetscInt));
1413: ISRestoreIndices(isrl,&ind);
1414: ISDestroy(&isrl);
1415: for (i=0; i<localsize; i++) *(indc+i) += size;
1416: ISCreateGeneral(PetscObjectComm((PetscObject)isy),localsize,indc,PETSC_OWN_POINTER,&isr);
1417: PetscObjectReference((PetscObject)isr);
1418: ISDestroy(&ilink->is);
1419: ilink->is = isr;
1420: PetscObjectReference((PetscObject)isr);
1421: ISDestroy(&ilink->is_col);
1422: ilink->is_col = isr;
1423: ISDestroy(&isr);
1424: KSPGetPC(ilink->ksp, &subpc);
1425: PetscObjectTypeCompare((PetscObject)subpc,PCFIELDSPLIT,&flg);
1426: if(flg) {
1427: IS iszl,isz;
1428: MPI_Comm comm;
1429: ISGetLocalSize(ilink->is,&localsize);
1430: comm = PetscObjectComm((PetscObject)ilink->is);
1431: ISEmbed(isy, ilink->is, PETSC_TRUE, &iszl);
1432: MPI_Scan(&localsize,&sizez,1,MPIU_INT,MPI_SUM,comm);
1433: sizez -= localsize;
1434: ISGetLocalSize(iszl,&localsize);
1435: PetscMalloc1(localsize,&indcz);
1436: ISGetIndices(iszl,&indz);
1437: PetscMemcpy(indcz,indz,localsize*sizeof(PetscInt));
1438: ISRestoreIndices(iszl,&indz);
1439: ISDestroy(&iszl);
1440: for (i=0; i<localsize; i++) *(indcz+i) += sizez;
1441: ISCreateGeneral(comm,localsize,indcz,PETSC_OWN_POINTER,&isz);
1442: PCFieldSplitRestrictIS(subpc,isz);
1443: ISDestroy(&isz);
1444: }
1445: next = ilink->next;
1446: ilink = next;
1447: }
1448: jac->isrestrict = PETSC_TRUE;
1449: return(0);
1450: }
1452: static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is)
1453: {
1454: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1455: PetscErrorCode ierr;
1456: PC_FieldSplitLink ilink, next = jac->head;
1457: char prefix[128];
1460: if (jac->splitdefined) {
1461: PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1462: return(0);
1463: }
1464: PetscNew(&ilink);
1465: if (splitname) {
1466: PetscStrallocpy(splitname,&ilink->splitname);
1467: } else {
1468: PetscMalloc1(8,&ilink->splitname);
1469: PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);
1470: }
1471: ilink->event = jac->nsplits < 5 ? KSP_Solve_FS_0 + jac->nsplits : KSP_Solve_FS_0 + 4; /* Any split great than 4 gets logged in the 4th split */
1472: PetscObjectReference((PetscObject)is);
1473: ISDestroy(&ilink->is);
1474: ilink->is = is;
1475: PetscObjectReference((PetscObject)is);
1476: ISDestroy(&ilink->is_col);
1477: ilink->is_col = is;
1478: ilink->next = NULL;
1479: KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1480: KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);
1481: PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1482: KSPSetType(ilink->ksp,KSPPREONLY);
1483: PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);
1485: PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);
1486: KSPSetOptionsPrefix(ilink->ksp,prefix);
1488: if (!next) {
1489: jac->head = ilink;
1490: ilink->previous = NULL;
1491: } else {
1492: while (next->next) {
1493: next = next->next;
1494: }
1495: next->next = ilink;
1496: ilink->previous = next;
1497: }
1498: jac->nsplits++;
1499: return(0);
1500: }
1502: /*@C
1503: PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner
1505: Logically Collective on PC
1507: Input Parameters:
1508: + pc - the preconditioner context
1509: . splitname - name of this split, if NULL the number of the split is used
1510: . n - the number of fields in this split
1511: - fields - the fields in this split
1513: Level: intermediate
1515: Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field.
1517: The PCFieldSplitSetFields() is for defining fields as strided blocks. For example, if the block
1518: size is three then one can define a field as 0, or 1 or 2 or 0,1 or 0,2 or 1,2 which mean
1519: 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x....
1520: where the numbered entries indicate what is in the field.
1522: This function is called once per split (it creates a new split each time). Solve options
1523: for this split will be available under the prefix -fieldsplit_SPLITNAME_.
1525: Developer Note: This routine does not actually create the IS representing the split, that is delayed
1526: until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be
1527: available when this routine is called.
1529: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS()
1531: @*/
1532: PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1533: {
1539: if (n < 1) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname);
1541: PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt*,const PetscInt*),(pc,splitname,n,fields,fields_col));
1542: return(0);
1543: }
1545: /*@
1546: PCFieldSplitSetDiagUseAmat - set flag indicating whether to extract diagonal blocks from Amat (rather than Pmat)
1548: Logically Collective on PC
1550: Input Parameters:
1551: + pc - the preconditioner object
1552: - flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from
1554: Options Database:
1555: . -pc_fieldsplit_diag_use_amat
1557: Level: intermediate
1559: .seealso: PCFieldSplitGetDiagUseAmat(), PCFieldSplitSetOffDiagUseAmat(), PCFIELDSPLIT
1561: @*/
1562: PetscErrorCode PCFieldSplitSetDiagUseAmat(PC pc,PetscBool flg)
1563: {
1564: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1565: PetscBool isfs;
1570: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1571: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1572: jac->diag_use_amat = flg;
1573: return(0);
1574: }
1576: /*@
1577: PCFieldSplitGetDiagUseAmat - get the flag indicating whether to extract diagonal blocks from Amat (rather than Pmat)
1579: Logically Collective on PC
1581: Input Parameters:
1582: . pc - the preconditioner object
1584: Output Parameters:
1585: . flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from
1588: Level: intermediate
1590: .seealso: PCFieldSplitSetDiagUseAmat(), PCFieldSplitGetOffDiagUseAmat(), PCFIELDSPLIT
1592: @*/
1593: PetscErrorCode PCFieldSplitGetDiagUseAmat(PC pc,PetscBool *flg)
1594: {
1595: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1596: PetscBool isfs;
1602: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1603: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1604: *flg = jac->diag_use_amat;
1605: return(0);
1606: }
1608: /*@
1609: PCFieldSplitSetOffDiagUseAmat - set flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat)
1611: Logically Collective on PC
1613: Input Parameters:
1614: + pc - the preconditioner object
1615: - flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from
1617: Options Database:
1618: . -pc_fieldsplit_off_diag_use_amat
1620: Level: intermediate
1622: .seealso: PCFieldSplitGetOffDiagUseAmat(), PCFieldSplitSetDiagUseAmat(), PCFIELDSPLIT
1624: @*/
1625: PetscErrorCode PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg)
1626: {
1627: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1628: PetscBool isfs;
1633: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1634: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1635: jac->offdiag_use_amat = flg;
1636: return(0);
1637: }
1639: /*@
1640: PCFieldSplitGetOffDiagUseAmat - get the flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat)
1642: Logically Collective on PC
1644: Input Parameters:
1645: . pc - the preconditioner object
1647: Output Parameters:
1648: . flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from
1651: Level: intermediate
1653: .seealso: PCFieldSplitSetOffDiagUseAmat(), PCFieldSplitGetDiagUseAmat(), PCFIELDSPLIT
1655: @*/
1656: PetscErrorCode PCFieldSplitGetOffDiagUseAmat(PC pc,PetscBool *flg)
1657: {
1658: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1659: PetscBool isfs;
1665: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1666: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1667: *flg = jac->offdiag_use_amat;
1668: return(0);
1669: }
1673: /*@C
1674: PCFieldSplitSetIS - Sets the exact elements for field
1676: Logically Collective on PC
1678: Input Parameters:
1679: + pc - the preconditioner context
1680: . splitname - name of this split, if NULL the number of the split is used
1681: - is - the index set that defines the vector elements in this field
1684: Notes:
1685: Use PCFieldSplitSetFields(), for fields defined by strided types.
1687: This function is called once per split (it creates a new split each time). Solve options
1688: for this split will be available under the prefix -fieldsplit_SPLITNAME_.
1690: Level: intermediate
1692: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize()
1694: @*/
1695: PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is)
1696: {
1703: PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));
1704: return(0);
1705: }
1707: /*@C
1708: PCFieldSplitGetIS - Retrieves the elements for a field as an IS
1710: Logically Collective on PC
1712: Input Parameters:
1713: + pc - the preconditioner context
1714: - splitname - name of this split
1716: Output Parameter:
1717: - is - the index set that defines the vector elements in this field, or NULL if the field is not found
1719: Level: intermediate
1721: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS()
1723: @*/
1724: PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is)
1725: {
1732: {
1733: PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;
1734: PC_FieldSplitLink ilink = jac->head;
1735: PetscBool found;
1737: *is = NULL;
1738: while (ilink) {
1739: PetscStrcmp(ilink->splitname, splitname, &found);
1740: if (found) {
1741: *is = ilink->is;
1742: break;
1743: }
1744: ilink = ilink->next;
1745: }
1746: }
1747: return(0);
1748: }
1750: /*@
1751: PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the
1752: fieldsplit preconditioner. If not set the matrix block size is used.
1754: Logically Collective on PC
1756: Input Parameters:
1757: + pc - the preconditioner context
1758: - bs - the block size
1760: Level: intermediate
1762: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields()
1764: @*/
1765: PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs)
1766: {
1772: PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));
1773: return(0);
1774: }
1776: /*@C
1777: PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits
1779: Collective on KSP
1781: Input Parameter:
1782: . pc - the preconditioner context
1784: Output Parameters:
1785: + n - the number of splits
1786: - subksp - the array of KSP contexts
1788: Note:
1789: After PCFieldSplitGetSubKSP() the array of KSPs is to be freed by the user with PetscFree()
1790: (not the KSP just the array that contains them).
1792: You must call KSPSetUp() before calling PCFieldSplitGetSubKSP().
1794: Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs.
1795: You can call PCFieldSplitGetSubKSP(pc,n,PETSC_NULL_KSP,ierr) to determine how large the
1796: KSP array must be.
1799: Level: advanced
1801: .seealso: PCFIELDSPLIT
1802: @*/
1803: PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
1804: {
1810: PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));
1811: return(0);
1812: }
1814: /*@
1815: PCFieldSplitSetSchurPre - Indicates what operator is used to construct the preconditioner for the Schur complement.
1816: A11 matrix. Otherwise no preconditioner is used.
1818: Collective on PC
1820: Input Parameters:
1821: + pc - the preconditioner context
1822: . ptype - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11 (default), PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_SCHUR_PRE_USER
1823: PC_FIELDSPLIT_SCHUR_PRE_SELFP, and PC_FIELDSPLIT_SCHUR_PRE_FULL
1824: - userpre - matrix to use for preconditioning, or NULL
1826: Options Database:
1827: . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11. See notes for meaning of various arguments
1829: Notes:
1830: $ If ptype is
1831: $ a11 then the preconditioner for the Schur complement is generated from the block diagonal part of the preconditioner
1832: $ matrix associated with the Schur complement (i.e. A11), not the Schur complement matrix
1833: $ self the preconditioner for the Schur complement is generated from the symbolic representation of the Schur complement matrix:
1834: $ The only preconditioner that currently works with this symbolic respresentation matrix object is the PCLSC
1835: $ preconditioner
1836: $ user then the preconditioner for the Schur complement is generated from the user provided matrix (pre argument
1837: $ to this function).
1838: $ selfp then the preconditioning for the Schur complement is generated from an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01
1839: $ This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be
1840: $ lumped before extracting the diagonal using the additional option -fieldsplit_1_mat_schur_complement_ainv_type lump
1841: $ full then the preconditioner for the Schur complement is generated from the exact Schur complement matrix representation computed internally by PFIELDSPLIT (this is expensive)
1842: $ useful mostly as a test that the Schur complement approach can work for your problem
1844: When solving a saddle point problem, where the A11 block is identically zero, using a11 as the ptype only makes sense
1845: with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and
1846: -fieldsplit_1_pc_type lsc which uses the least squares commutator to compute a preconditioner for the Schur complement.
1848: Level: intermediate
1850: .seealso: PCFieldSplitGetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType,
1851: MatSchurComplementSetAinvType(), PCLSC
1853: @*/
1854: PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
1855: {
1860: PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));
1861: return(0);
1862: }
1863: PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) {return PCFieldSplitSetSchurPre(pc,ptype,pre);} /* Deprecated name */
1865: /*@
1866: PCFieldSplitGetSchurPre - For Schur complement fieldsplit, determine how the Schur complement will be
1867: preconditioned. See PCFieldSplitSetSchurPre() for details.
1869: Logically Collective on PC
1871: Input Parameters:
1872: . pc - the preconditioner context
1874: Output Parameters:
1875: + ptype - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11, PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_PRE_USER
1876: - userpre - matrix to use for preconditioning (with PC_FIELDSPLIT_PRE_USER), or NULL
1878: Level: intermediate
1880: .seealso: PCFieldSplitSetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC
1882: @*/
1883: PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
1884: {
1889: PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre));
1890: return(0);
1891: }
1893: /*@
1894: PCFieldSplitSchurGetS - extract the MatSchurComplement object used by this PC in case it needs to be configured separately
1896: Not collective
1898: Input Parameter:
1899: . pc - the preconditioner context
1901: Output Parameter:
1902: . S - the Schur complement matrix
1904: Notes:
1905: This matrix should not be destroyed using MatDestroy(); rather, use PCFieldSplitSchurRestoreS().
1907: Level: advanced
1909: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurRestoreS()
1911: @*/
1912: PetscErrorCode PCFieldSplitSchurGetS(PC pc,Mat *S)
1913: {
1915: const char* t;
1916: PetscBool isfs;
1917: PC_FieldSplit *jac;
1921: PetscObjectGetType((PetscObject)pc,&t);
1922: PetscStrcmp(t,PCFIELDSPLIT,&isfs);
1923: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
1924: jac = (PC_FieldSplit*)pc->data;
1925: if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
1926: if (S) *S = jac->schur;
1927: return(0);
1928: }
1930: /*@
1931: PCFieldSplitSchurRestoreS - restores the MatSchurComplement object used by this PC
1933: Not collective
1935: Input Parameters:
1936: + pc - the preconditioner context
1937: . S - the Schur complement matrix
1939: Level: advanced
1941: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurGetS()
1943: @*/
1944: PetscErrorCode PCFieldSplitSchurRestoreS(PC pc,Mat *S)
1945: {
1947: const char* t;
1948: PetscBool isfs;
1949: PC_FieldSplit *jac;
1953: PetscObjectGetType((PetscObject)pc,&t);
1954: PetscStrcmp(t,PCFIELDSPLIT,&isfs);
1955: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
1956: jac = (PC_FieldSplit*)pc->data;
1957: if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
1958: if (!S || *S != jac->schur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten");
1959: return(0);
1960: }
1963: static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
1964: {
1965: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1969: jac->schurpre = ptype;
1970: if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) {
1971: MatDestroy(&jac->schur_user);
1972: jac->schur_user = pre;
1973: PetscObjectReference((PetscObject)jac->schur_user);
1974: }
1975: return(0);
1976: }
1978: static PetscErrorCode PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
1979: {
1980: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1983: *ptype = jac->schurpre;
1984: *pre = jac->schur_user;
1985: return(0);
1986: }
1988: /*@
1989: PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain in the preconditioner
1991: Collective on PC
1993: Input Parameters:
1994: + pc - the preconditioner context
1995: - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default
1997: Options Database:
1998: . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full
2001: Level: intermediate
2003: Notes:
2004: The FULL factorization is
2006: $ (A B) = (1 0) (A 0) (1 Ainv*B) = L D U
2007: $ (C E) (C*Ainv 1) (0 S) (0 1 )
2009: where S = E - C*Ainv*B. In practice, the full factorization is applied via block triangular solves with the grouping L*(D*U). UPPER uses D*U, LOWER uses L*D,
2010: and DIAG is the diagonal part with the sign of S flipped (because this makes the preconditioner positive definite for many formulations, thus allowing the use of KSPMINRES). Sign flipping of S can be turned off with PCFieldSplitSetSchurScale().
2012: $ If A and S are solved exactly
2013: $ *) FULL factorization is a direct solver.
2014: $ *) The preconditioned operator with LOWER or UPPER has all eigenvalues equal to 1 and minimal polynomial of degree 2, so KSPGMRES converges in 2 iterations.
2015: $ *) With DIAG, the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so KSPGMRES converges in at most 4 iterations.
2017: If the iteration count is very low, consider using KSPFGMRES or KSPGCR which can use one less preconditioner
2018: application in this case. Note that the preconditioned operator may be highly non-normal, so such fast convergence may not be observed in practice.
2020: For symmetric problems in which A is positive definite and S is negative definite, DIAG can be used with KSPMINRES.
2022: Note that a flexible method like KSPFGMRES or KSPGCR must be used if the fieldsplit preconditioner is nonlinear (e.g. a few iterations of a Krylov method is used to solve with A or S).
2024: References:
2025: + 1. - Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000).
2026: - 2. - Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001).
2028: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCFieldSplitSetSchurScale()
2029: @*/
2030: PetscErrorCode PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype)
2031: {
2036: PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));
2037: return(0);
2038: }
2040: static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype)
2041: {
2042: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2045: jac->schurfactorization = ftype;
2046: return(0);
2047: }
2049: /*@
2050: PCFieldSplitSetSchurScale - Controls the sign flip of S for PC_FIELDSPLIT_SCHUR_FACT_DIAG.
2052: Collective on PC
2054: Input Parameters:
2055: + pc - the preconditioner context
2056: - scale - scaling factor for the Schur complement
2058: Options Database:
2059: . -pc_fieldsplit_schur_scale - default is -1.0
2061: Level: intermediate
2063: .seealso: PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurFactType, PCFieldSplitSetSchurScale()
2064: @*/
2065: PetscErrorCode PCFieldSplitSetSchurScale(PC pc,PetscScalar scale)
2066: {
2072: PetscTryMethod(pc,"PCFieldSplitSetSchurScale_C",(PC,PetscScalar),(pc,scale));
2073: return(0);
2074: }
2076: static PetscErrorCode PCFieldSplitSetSchurScale_FieldSplit(PC pc,PetscScalar scale)
2077: {
2078: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2081: jac->schurscale = scale;
2082: return(0);
2083: }
2085: /*@C
2086: PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement
2088: Collective on KSP
2090: Input Parameter:
2091: . pc - the preconditioner context
2093: Output Parameters:
2094: + A00 - the (0,0) block
2095: . A01 - the (0,1) block
2096: . A10 - the (1,0) block
2097: - A11 - the (1,1) block
2099: Level: advanced
2101: .seealso: PCFIELDSPLIT
2102: @*/
2103: PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11)
2104: {
2105: PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;
2109: if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach.");
2110: if (A00) *A00 = jac->pmat[0];
2111: if (A01) *A01 = jac->B;
2112: if (A10) *A10 = jac->C;
2113: if (A11) *A11 = jac->pmat[1];
2114: return(0);
2115: }
2117: static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type)
2118: {
2119: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2123: jac->type = type;
2124: if (type == PC_COMPOSITE_SCHUR) {
2125: pc->ops->apply = PCApply_FieldSplit_Schur;
2126: pc->ops->view = PCView_FieldSplit_Schur;
2128: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);
2129: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit);
2130: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit);
2131: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);
2132: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",PCFieldSplitSetSchurScale_FieldSplit);
2134: } else {
2135: pc->ops->apply = PCApply_FieldSplit;
2136: pc->ops->view = PCView_FieldSplit;
2138: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2139: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",0);
2140: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",0);
2141: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);
2142: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",0);
2143: }
2144: return(0);
2145: }
2147: static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs)
2148: {
2149: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2152: if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs);
2153: if (jac->bs > 0 && jac->bs != bs) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Cannot change fieldsplit blocksize from %D to %D after it has been set",jac->bs,bs);
2154: jac->bs = bs;
2155: return(0);
2156: }
2158: /*@
2159: PCFieldSplitSetType - Sets the type of fieldsplit preconditioner.
2161: Collective on PC
2163: Input Parameter:
2164: . pc - the preconditioner context
2165: . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR
2167: Options Database Key:
2168: . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type
2170: Level: Intermediate
2172: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
2174: .seealso: PCCompositeSetType()
2176: @*/
2177: PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type)
2178: {
2183: PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));
2184: return(0);
2185: }
2187: /*@
2188: PCFieldSplitGetType - Gets the type of fieldsplit preconditioner.
2190: Not collective
2192: Input Parameter:
2193: . pc - the preconditioner context
2195: Output Parameter:
2196: . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR
2198: Level: Intermediate
2200: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
2201: .seealso: PCCompositeSetType()
2202: @*/
2203: PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type)
2204: {
2205: PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;
2210: *type = jac->type;
2211: return(0);
2212: }
2214: /*@
2215: PCFieldSplitSetDMSplits - Flags whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible.
2217: Logically Collective
2219: Input Parameters:
2220: + pc - the preconditioner context
2221: - flg - boolean indicating whether to use field splits defined by the DM
2223: Options Database Key:
2224: . -pc_fieldsplit_dm_splits
2226: Level: Intermediate
2228: .keywords: PC, DM, composite preconditioner, additive, multiplicative
2230: .seealso: PCFieldSplitGetDMSplits()
2232: @*/
2233: PetscErrorCode PCFieldSplitSetDMSplits(PC pc,PetscBool flg)
2234: {
2235: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2236: PetscBool isfs;
2242: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2243: if (isfs) {
2244: jac->dm_splits = flg;
2245: }
2246: return(0);
2247: }
2250: /*@
2251: PCFieldSplitGetDMSplits - Returns flag indicating whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible.
2253: Logically Collective
2255: Input Parameter:
2256: . pc - the preconditioner context
2258: Output Parameter:
2259: . flg - boolean indicating whether to use field splits defined by the DM
2261: Level: Intermediate
2263: .keywords: PC, DM, composite preconditioner, additive, multiplicative
2265: .seealso: PCFieldSplitSetDMSplits()
2267: @*/
2268: PetscErrorCode PCFieldSplitGetDMSplits(PC pc,PetscBool* flg)
2269: {
2270: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2271: PetscBool isfs;
2277: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2278: if (isfs) {
2279: if(flg) *flg = jac->dm_splits;
2280: }
2281: return(0);
2282: }
2284: /* -------------------------------------------------------------------------------------*/
2285: /*MC
2286: PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual
2287: fields or groups of fields. See the users manual section "Solving Block Matrices" for more details.
2289: To set options on the solvers for each block append -fieldsplit_ to all the PC
2290: options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1
2292: To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP()
2293: and set the options directly on the resulting KSP object
2295: Level: intermediate
2297: Options Database Keys:
2298: + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split
2299: . -pc_fieldsplit_default - automatically add any fields to additional splits that have not
2300: been supplied explicitly by -pc_fieldsplit_%d_fields
2301: . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields)
2302: . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur> - type of relaxation or factorization splitting
2303: . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11; see PCFieldSplitSetSchurPre()
2304: . -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero or negative diagonal and uses Schur complement with no preconditioner as the solver
2306: - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_
2307: for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields
2309: Notes:
2310: Use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS()
2311: to define a field by an arbitrary collection of entries.
2313: If no fields are set the default is used. The fields are defined by entries strided by bs,
2314: beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(),
2315: if this is not called the block size defaults to the blocksize of the second matrix passed
2316: to KSPSetOperators()/PCSetOperators().
2318: $ For the Schur complement preconditioner if J = ( A00 A01 )
2319: $ ( A10 A11 )
2320: $ the preconditioner using full factorization is
2321: $ ( I -ksp(A00) A01 ) ( inv(A00) 0 ) ( I 0 )
2322: $ ( 0 I ) ( 0 ksp(S) ) ( -A10 ksp(A00) I )
2323: where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. S is the Schur complement
2324: $ S = A11 - A10 ksp(A00) A01
2325: which is usually dense and not stored explicitly. The action of ksp(S) is computed using the KSP solver with prefix -fieldsplit_splitname_ (where splitname was given
2326: in providing the SECOND split or 1 if not give). For PCFieldSplitGetKSP() when field number is 0,
2327: it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default
2328: A11 is used to construct a preconditioner for S, use PCFieldSplitSetSchurPre() for all the possible ways to construct the preconditioner for S.
2330: The factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above,
2331: diag gives
2332: $ ( inv(A00) 0 )
2333: $ ( 0 -ksp(S) )
2334: note that slightly counter intuitively there is a negative in front of the ksp(S) so that the preconditioner is positive definite. For SPD matrices J, the sign flip
2335: can be turned off with PCFieldSplitSetSchurScale() or by command line -pc_fieldsplit_schur_scale 1.0. The lower factorization is the inverse of
2336: $ ( A00 0 )
2337: $ ( A10 S )
2338: where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of
2339: $ ( A00 A01 )
2340: $ ( 0 S )
2341: where again the inverses of A00 and S are applied using KSPs.
2343: If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS
2344: is used automatically for a second block.
2346: The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1.
2347: Generally it should be used with the AIJ format.
2349: The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see,
2350: for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT
2351: inside a smoother resulting in "Distributive Smoothers".
2353: Concepts: physics based preconditioners, block preconditioners
2355: There is a nice discussion of block preconditioners in
2357: [El08] A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier-Stokes equations
2358: Howard Elman, V.E. Howle, John Shadid, Robert Shuttleworth, Ray Tuminaro, Journal of Computational Physics 227 (2008) 1790--1808
2359: http://chess.cs.umd.edu/~elman/papers/tax.pdf
2361: The Constrained Pressure Preconditioner (CPR) can be implemented using PCCOMPOSITE with PCGALERKIN. CPR first solves an R A P subsystem, updates the
2362: residual on all variables (PCCompositeSetType(pc,PC_COMPOSITE_MULTIPLICATIVE)), and then applies a simple ILU like preconditioner on all the variables.
2364: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC,
2365: PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre(),
2366: MatSchurComplementSetAinvType(), PCFieldSplitSetSchurScale()
2367: M*/
2369: PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc)
2370: {
2372: PC_FieldSplit *jac;
2375: PetscNewLog(pc,&jac);
2377: jac->bs = -1;
2378: jac->nsplits = 0;
2379: jac->type = PC_COMPOSITE_MULTIPLICATIVE;
2380: jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */
2381: jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL;
2382: jac->schurscale = -1.0;
2383: jac->dm_splits = PETSC_TRUE;
2385: pc->data = (void*)jac;
2387: pc->ops->apply = PCApply_FieldSplit;
2388: pc->ops->applytranspose = PCApplyTranspose_FieldSplit;
2389: pc->ops->setup = PCSetUp_FieldSplit;
2390: pc->ops->reset = PCReset_FieldSplit;
2391: pc->ops->destroy = PCDestroy_FieldSplit;
2392: pc->ops->setfromoptions = PCSetFromOptions_FieldSplit;
2393: pc->ops->view = PCView_FieldSplit;
2394: pc->ops->applyrichardson = 0;
2396: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2397: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);
2398: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);
2399: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);
2400: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);
2401: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",PCFieldSplitRestrictIS_FieldSplit);
2402: return(0);
2403: }