Actual source code: fieldsplit.c
petsc-3.7.7 2017-09-25
2: #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/
3: #include <petsc/private/kspimpl.h>
4: #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,is_orig;
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: PC_FieldSplitLink head;
48: PetscBool reset; /* indicates PCReset() has been last called on this object, hack */
49: PetscBool isrestrict; /* indicates PCFieldSplitRestrictIS() has been last called on this object, hack */
50: PetscBool suboptionsset; /* Indicates that the KSPSetFromOptions() has been called on the sub-KSPs */
51: PetscBool dm_splits; /* Whether to use DMCreateFieldDecomposition() whenever possible */
52: PetscBool diag_use_amat; /* Whether to extract diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */
53: PetscBool offdiag_use_amat; /* Whether to extract off-diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */
54: } PC_FieldSplit;
56: /*
57: Notes: there is no particular reason that pmat, x, and y are stored as arrays in PC_FieldSplit instead of
58: inside PC_FieldSplitLink, just historical. If you want to be able to add new fields after already using the
59: PC you could change this.
60: */
62: /* This helper is so that setting a user-provided preconditioning matrix is orthogonal to choosing to use it. This way the
63: * application-provided FormJacobian can provide this matrix without interfering with the user's (command-line) choices. */
64: static Mat FieldSplitSchurPre(PC_FieldSplit *jac)
65: {
66: switch (jac->schurpre) {
67: case PC_FIELDSPLIT_SCHUR_PRE_SELF: return jac->schur;
68: case PC_FIELDSPLIT_SCHUR_PRE_SELFP: return jac->schurp;
69: case PC_FIELDSPLIT_SCHUR_PRE_A11: return jac->pmat[1];
70: case PC_FIELDSPLIT_SCHUR_PRE_FULL: /* We calculate this and store it in schur_user */
71: case PC_FIELDSPLIT_SCHUR_PRE_USER: /* Use a user-provided matrix if it is given, otherwise diagonal block */
72: default:
73: return jac->schur_user ? jac->schur_user : jac->pmat[1];
74: }
75: }
78: #include <petscdraw.h>
81: static PetscErrorCode PCView_FieldSplit(PC pc,PetscViewer viewer)
82: {
83: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
84: PetscErrorCode ierr;
85: PetscBool iascii,isdraw;
86: PetscInt i,j;
87: PC_FieldSplitLink ilink = jac->head;
90: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
91: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
92: if (iascii) {
93: if (jac->bs > 0) {
94: PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs);
95: } else {
96: PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D\n",PCCompositeTypes[jac->type],jac->nsplits);
97: }
98: if (pc->useAmat) {
99: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");
100: }
101: if (jac->diag_use_amat) {
102: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for diagonal blocks\n");
103: }
104: if (jac->offdiag_use_amat) {
105: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for off-diagonal blocks\n");
106: }
107: PetscViewerASCIIPrintf(viewer," Solver info for each split is in the following KSP objects:\n");
108: PetscViewerASCIIPushTab(viewer);
109: for (i=0; i<jac->nsplits; i++) {
110: if (ilink->fields) {
111: PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);
112: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
113: for (j=0; j<ilink->nfields; j++) {
114: if (j > 0) {
115: PetscViewerASCIIPrintf(viewer,",");
116: }
117: PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);
118: }
119: PetscViewerASCIIPrintf(viewer,"\n");
120: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
121: } else {
122: PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);
123: }
124: KSPView(ilink->ksp,viewer);
125: ilink = ilink->next;
126: }
127: PetscViewerASCIIPopTab(viewer);
128: }
130: if (isdraw) {
131: PetscDraw draw;
132: PetscReal x,y,w,wd;
134: PetscViewerDrawGetDraw(viewer,0,&draw);
135: PetscDrawGetCurrentPoint(draw,&x,&y);
136: w = 2*PetscMin(1.0 - x,x);
137: wd = w/(jac->nsplits + 1);
138: x = x - wd*(jac->nsplits-1)/2.0;
139: for (i=0; i<jac->nsplits; i++) {
140: PetscDrawPushCurrentPoint(draw,x,y);
141: KSPView(ilink->ksp,viewer);
142: PetscDrawPopCurrentPoint(draw);
143: x += wd;
144: ilink = ilink->next;
145: }
146: }
147: return(0);
148: }
152: static PetscErrorCode PCView_FieldSplit_Schur(PC pc,PetscViewer viewer)
153: {
154: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
155: PetscErrorCode ierr;
156: PetscBool iascii,isdraw;
157: PetscInt i,j;
158: PC_FieldSplitLink ilink = jac->head;
161: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
162: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
163: if (iascii) {
164: if (jac->bs > 0) {
165: PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, blocksize = %D, factorization %s\n",jac->bs,PCFieldSplitSchurFactTypes[jac->schurfactorization]);
166: } else {
167: PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, factorization %s\n",PCFieldSplitSchurFactTypes[jac->schurfactorization]);
168: }
169: if (pc->useAmat) {
170: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");
171: }
172: switch (jac->schurpre) {
173: case PC_FIELDSPLIT_SCHUR_PRE_SELF:
174: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from S itself\n");break;
175: case PC_FIELDSPLIT_SCHUR_PRE_SELFP:
176: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from Sp, an assembled approximation to S, which uses (lumped, if requested) A00's diagonal's inverse\n");break;
177: case PC_FIELDSPLIT_SCHUR_PRE_A11:
178: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n");break;
179: case PC_FIELDSPLIT_SCHUR_PRE_FULL:
180: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from the exact Schur complement\n");break;
181: case PC_FIELDSPLIT_SCHUR_PRE_USER:
182: if (jac->schur_user) {
183: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from user provided matrix\n");
184: } else {
185: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n");
186: }
187: break;
188: default:
189: SETERRQ1(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Invalid Schur preconditioning type: %d", jac->schurpre);
190: }
191: PetscViewerASCIIPrintf(viewer," Split info:\n");
192: PetscViewerASCIIPushTab(viewer);
193: for (i=0; i<jac->nsplits; i++) {
194: if (ilink->fields) {
195: PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);
196: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
197: for (j=0; j<ilink->nfields; j++) {
198: if (j > 0) {
199: PetscViewerASCIIPrintf(viewer,",");
200: }
201: PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);
202: }
203: PetscViewerASCIIPrintf(viewer,"\n");
204: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
205: } else {
206: PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);
207: }
208: ilink = ilink->next;
209: }
210: PetscViewerASCIIPrintf(viewer,"KSP solver for A00 block\n");
211: PetscViewerASCIIPushTab(viewer);
212: if (jac->head) {
213: KSPView(jac->head->ksp,viewer);
214: } else {PetscViewerASCIIPrintf(viewer," not yet available\n");}
215: PetscViewerASCIIPopTab(viewer);
216: if (jac->head && jac->kspupper != jac->head->ksp) {
217: PetscViewerASCIIPrintf(viewer,"KSP solver for upper A00 in upper triangular factor \n");
218: PetscViewerASCIIPushTab(viewer);
219: if (jac->kspupper) {KSPView(jac->kspupper,viewer);}
220: else {PetscViewerASCIIPrintf(viewer," not yet available\n");}
221: PetscViewerASCIIPopTab(viewer);
222: }
223: PetscViewerASCIIPrintf(viewer,"KSP solver for S = A11 - A10 inv(A00) A01 \n");
224: PetscViewerASCIIPushTab(viewer);
225: if (jac->kspschur) {
226: KSPView(jac->kspschur,viewer);
227: } else {
228: PetscViewerASCIIPrintf(viewer," not yet available\n");
229: }
230: PetscViewerASCIIPopTab(viewer);
231: PetscViewerASCIIPopTab(viewer);
232: } else if (isdraw && jac->head) {
233: PetscDraw draw;
234: PetscReal x,y,w,wd,h;
235: PetscInt cnt = 2;
236: char str[32];
238: PetscViewerDrawGetDraw(viewer,0,&draw);
239: PetscDrawGetCurrentPoint(draw,&x,&y);
240: if (jac->kspupper != jac->head->ksp) cnt++;
241: w = 2*PetscMin(1.0 - x,x);
242: wd = w/(cnt + 1);
244: PetscSNPrintf(str,32,"Schur fact. %s",PCFieldSplitSchurFactTypes[jac->schurfactorization]);
245: PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
246: y -= h;
247: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_USER && !jac->schur_user) {
248: PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[PC_FIELDSPLIT_SCHUR_PRE_A11]);
249: } else {
250: PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[jac->schurpre]);
251: }
252: PetscDrawStringBoxed(draw,x+wd*(cnt-1)/2.0,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
253: y -= h;
254: x = x - wd*(cnt-1)/2.0;
256: PetscDrawPushCurrentPoint(draw,x,y);
257: KSPView(jac->head->ksp,viewer);
258: PetscDrawPopCurrentPoint(draw);
259: if (jac->kspupper != jac->head->ksp) {
260: x += wd;
261: PetscDrawPushCurrentPoint(draw,x,y);
262: KSPView(jac->kspupper,viewer);
263: PetscDrawPopCurrentPoint(draw);
264: }
265: x += wd;
266: PetscDrawPushCurrentPoint(draw,x,y);
267: KSPView(jac->kspschur,viewer);
268: PetscDrawPopCurrentPoint(draw);
269: }
270: return(0);
271: }
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: }
318: static PetscErrorCode PCFieldSplitSetDefaults(PC pc)
319: {
320: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
321: PetscErrorCode ierr;
322: PC_FieldSplitLink ilink = jac->head;
323: PetscBool fieldsplit_default = PETSC_FALSE,stokes = PETSC_FALSE,coupling = PETSC_FALSE;
324: PetscInt i;
327: /*
328: Kinda messy, but at least this now uses DMCreateFieldDecomposition() even with jac->reset.
329: Should probably be rewritten.
330: */
331: if (!ilink || jac->reset) {
332: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);
333: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_coupling",&coupling,NULL);
334: if (pc->dm && jac->dm_splits && !stokes && !coupling) {
335: PetscInt numFields, f, i, j;
336: char **fieldNames;
337: IS *fields;
338: DM *dms;
339: DM subdm[128];
340: PetscBool flg;
342: DMCreateFieldDecomposition(pc->dm, &numFields, &fieldNames, &fields, &dms);
343: /* Allow the user to prescribe the splits */
344: for (i = 0, flg = PETSC_TRUE;; i++) {
345: PetscInt ifields[128];
346: IS compField;
347: char optionname[128], splitname[8];
348: PetscInt nfields = numFields;
350: PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%D_fields", i);
351: PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix, optionname, ifields, &nfields, &flg);
352: if (!flg) break;
353: if (numFields > 128) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot currently support %d > 128 fields", numFields);
354: DMCreateSubDM(pc->dm, nfields, ifields, &compField, &subdm[i]);
355: if (nfields == 1) {
356: PCFieldSplitSetIS(pc, fieldNames[ifields[0]], compField);
357: /* PetscPrintf(PetscObjectComm((PetscObject)pc), "%s Field Indices:", fieldNames[ifields[0]]);
358: ISView(compField, NULL); */
359: } else {
360: PetscSNPrintf(splitname, sizeof(splitname), "%D", i);
361: PCFieldSplitSetIS(pc, splitname, compField);
362: /* PetscPrintf(PetscObjectComm((PetscObject)pc), "%s Field Indices:", splitname);
363: ISView(compField, NULL); */
364: }
365: ISDestroy(&compField);
366: for (j = 0; j < nfields; ++j) {
367: f = ifields[j];
368: PetscFree(fieldNames[f]);
369: ISDestroy(&fields[f]);
370: }
371: }
372: if (i == 0) {
373: for (f = 0; f < numFields; ++f) {
374: PCFieldSplitSetIS(pc, fieldNames[f], fields[f]);
375: PetscFree(fieldNames[f]);
376: ISDestroy(&fields[f]);
377: }
378: } else {
379: for (j=0; j<numFields; j++) {
380: DMDestroy(dms+j);
381: }
382: PetscFree(dms);
383: PetscMalloc1(i, &dms);
384: for (j = 0; j < i; ++j) dms[j] = subdm[j];
385: }
386: PetscFree(fieldNames);
387: PetscFree(fields);
388: if (dms) {
389: PetscInfo(pc, "Setting up physics based fieldsplit preconditioner using the embedded DM\n");
390: for (ilink = jac->head, i = 0; ilink; ilink = ilink->next, ++i) {
391: const char *prefix;
392: PetscObjectGetOptionsPrefix((PetscObject)(ilink->ksp),&prefix);
393: PetscObjectSetOptionsPrefix((PetscObject)(dms[i]), prefix);
394: KSPSetDM(ilink->ksp, dms[i]);
395: KSPSetDMActive(ilink->ksp, PETSC_FALSE);
396: PetscObjectIncrementTabLevel((PetscObject)dms[i],(PetscObject)ilink->ksp,0);
397: DMDestroy(&dms[i]);
398: }
399: PetscFree(dms);
400: }
401: } else {
402: if (jac->bs <= 0) {
403: if (pc->pmat) {
404: MatGetBlockSize(pc->pmat,&jac->bs);
405: } else jac->bs = 1;
406: }
408: if (stokes) {
409: IS zerodiags,rest;
410: PetscInt nmin,nmax;
412: MatGetOwnershipRange(pc->mat,&nmin,&nmax);
413: MatFindZeroDiagonals(pc->mat,&zerodiags);
414: ISComplement(zerodiags,nmin,nmax,&rest);
415: if (jac->reset) {
416: jac->head->is = rest;
417: jac->head->next->is = zerodiags;
418: } else {
419: PCFieldSplitSetIS(pc,"0",rest);
420: PCFieldSplitSetIS(pc,"1",zerodiags);
421: }
422: ISDestroy(&zerodiags);
423: ISDestroy(&rest);
424: } else if (coupling) {
425: IS coupling,rest;
426: PetscInt nmin,nmax;
428: MatGetOwnershipRange(pc->mat,&nmin,&nmax);
429: MatFindOffBlockDiagonalEntries(pc->mat,&coupling);
430: ISCreateStride(PetscObjectComm((PetscObject)pc->mat),nmax-nmin,nmin,1,&rest);
431: ISSetIdentity(rest);
432: if (jac->reset) {
433: jac->head->is = rest;
434: jac->head->next->is = coupling;
435: } else {
436: PCFieldSplitSetIS(pc,"0",rest);
437: PCFieldSplitSetIS(pc,"1",coupling);
438: }
439: ISDestroy(&coupling);
440: ISDestroy(&rest);
441: } else {
442: if (jac->reset && !jac->isrestrict) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cases not yet handled when PCReset() was used");
443: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&fieldsplit_default,NULL);
444: if (!fieldsplit_default) {
445: /* Allow user to set fields from command line, if bs was known at the time of PCSetFromOptions_FieldSplit()
446: then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */
447: PCFieldSplitSetRuntimeSplits_Private(pc);
448: if (jac->splitdefined) {PetscInfo(pc,"Splits defined using the options database\n");}
449: }
450: if ((fieldsplit_default || !jac->splitdefined) && !jac->isrestrict) {
451: PetscInfo(pc,"Using default splitting of fields\n");
452: for (i=0; i<jac->bs; i++) {
453: char splitname[8];
454: PetscSNPrintf(splitname,sizeof(splitname),"%D",i);
455: PCFieldSplitSetFields(pc,splitname,1,&i,&i);
456: }
457: jac->defaultsplit = PETSC_TRUE;
458: }
459: }
460: }
461: } else if (jac->nsplits == 1) {
462: if (ilink->is) {
463: IS is2;
464: PetscInt nmin,nmax;
466: MatGetOwnershipRange(pc->mat,&nmin,&nmax);
467: ISComplement(ilink->is,nmin,nmax,&is2);
468: PCFieldSplitSetIS(pc,"1",is2);
469: ISDestroy(&is2);
470: } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Must provide at least two sets of fields to PCFieldSplit()");
471: }
474: if (jac->nsplits < 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Unhandled case, must have at least two fields, not %d", jac->nsplits);
475: return(0);
476: }
478: PETSC_EXTERN PetscErrorCode PetscOptionsFindPairPrefix_Private(PetscOptions,const char pre[], const char name[], char *value[], PetscBool *flg);
482: static PetscErrorCode PCSetUp_FieldSplit(PC pc)
483: {
484: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
485: PetscErrorCode ierr;
486: PC_FieldSplitLink ilink;
487: PetscInt i,nsplit;
488: PetscBool sorted, sorted_col;
491: PCFieldSplitSetDefaults(pc);
492: nsplit = jac->nsplits;
493: ilink = jac->head;
495: /* get the matrices for each split */
496: if (!jac->issetup) {
497: PetscInt rstart,rend,nslots,bs;
499: jac->issetup = PETSC_TRUE;
501: /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */
502: if (jac->defaultsplit || !ilink->is) {
503: if (jac->bs <= 0) jac->bs = nsplit;
504: }
505: bs = jac->bs;
506: MatGetOwnershipRange(pc->pmat,&rstart,&rend);
507: nslots = (rend - rstart)/bs;
508: for (i=0; i<nsplit; i++) {
509: if (jac->defaultsplit) {
510: ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+i,nsplit,&ilink->is);
511: ISDuplicate(ilink->is,&ilink->is_col);
512: } else if (!ilink->is) {
513: if (ilink->nfields > 1) {
514: PetscInt *ii,*jj,j,k,nfields = ilink->nfields,*fields = ilink->fields,*fields_col = ilink->fields_col;
515: PetscMalloc1(ilink->nfields*nslots,&ii);
516: PetscMalloc1(ilink->nfields*nslots,&jj);
517: for (j=0; j<nslots; j++) {
518: for (k=0; k<nfields; k++) {
519: ii[nfields*j + k] = rstart + bs*j + fields[k];
520: jj[nfields*j + k] = rstart + bs*j + fields_col[k];
521: }
522: }
523: ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,ii,PETSC_OWN_POINTER,&ilink->is);
524: ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,jj,PETSC_OWN_POINTER,&ilink->is_col);
525: ISSetBlockSize(ilink->is, nfields);
526: ISSetBlockSize(ilink->is_col, nfields);
527: } else {
528: ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields[0],bs,&ilink->is);
529: ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields_col[0],bs,&ilink->is_col);
530: }
531: }
532: ISSorted(ilink->is,&sorted);
533: if (ilink->is_col) { ISSorted(ilink->is_col,&sorted_col); }
534: if (!sorted || !sorted_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Fields must be sorted when creating split");
535: ilink = ilink->next;
536: }
537: }
539: ilink = jac->head;
540: if (!jac->pmat) {
541: Vec xtmp;
543: MatCreateVecs(pc->pmat,&xtmp,NULL);
544: PetscMalloc1(nsplit,&jac->pmat);
545: PetscMalloc2(nsplit,&jac->x,nsplit,&jac->y);
546: for (i=0; i<nsplit; i++) {
547: MatNullSpace sp;
549: /* Check for preconditioning matrix attached to IS */
550: PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &jac->pmat[i]);
551: if (jac->pmat[i]) {
552: PetscObjectReference((PetscObject) jac->pmat[i]);
553: if (jac->type == PC_COMPOSITE_SCHUR) {
554: jac->schur_user = jac->pmat[i];
556: PetscObjectReference((PetscObject) jac->schur_user);
557: }
558: } else {
559: const char *prefix;
560: MatGetSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->pmat[i]);
561: KSPGetOptionsPrefix(ilink->ksp,&prefix);
562: MatSetOptionsPrefix(jac->pmat[i],prefix);
563: MatViewFromOptions(jac->pmat[i],NULL,"-mat_view");
564: }
565: /* create work vectors for each split */
566: MatCreateVecs(jac->pmat[i],&jac->x[i],&jac->y[i]);
567: ilink->x = jac->x[i]; ilink->y = jac->y[i]; ilink->z = NULL;
568: /* compute scatter contexts needed by multiplicative versions and non-default splits */
569: VecScatterCreate(xtmp,ilink->is,jac->x[i],NULL,&ilink->sctx);
570: PetscObjectQuery((PetscObject) ilink->is, "nearnullspace", (PetscObject*) &sp);
571: if (sp) {
572: MatSetNearNullSpace(jac->pmat[i], sp);
573: }
574: ilink = ilink->next;
575: }
576: VecDestroy(&xtmp);
577: } else {
578: for (i=0; i<nsplit; i++) {
579: Mat pmat;
581: /* Check for preconditioning matrix attached to IS */
582: PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &pmat);
583: if (!pmat) {
584: MatGetSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_REUSE_MATRIX,&jac->pmat[i]);
585: }
586: ilink = ilink->next;
587: }
588: }
589: if (jac->diag_use_amat) {
590: ilink = jac->head;
591: if (!jac->mat) {
592: PetscMalloc1(nsplit,&jac->mat);
593: for (i=0; i<nsplit; i++) {
594: MatGetSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->mat[i]);
595: ilink = ilink->next;
596: }
597: } else {
598: for (i=0; i<nsplit; i++) {
599: if (jac->mat[i]) {MatGetSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_REUSE_MATRIX,&jac->mat[i]);}
600: ilink = ilink->next;
601: }
602: }
603: } else {
604: jac->mat = jac->pmat;
605: }
607: /* Check for null space attached to IS */
608: ilink = jac->head;
609: for (i=0; i<nsplit; i++) {
610: MatNullSpace sp;
612: PetscObjectQuery((PetscObject) ilink->is, "nullspace", (PetscObject*) &sp);
613: if (sp) {
614: MatSetNullSpace(jac->mat[i], sp);
615: }
616: ilink = ilink->next;
617: }
619: if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR) {
620: /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */
621: /* FIXME: Can/should we reuse jac->mat whenever (jac->diag_use_amat) is true? */
622: ilink = jac->head;
623: if (nsplit == 2 && jac->type == PC_COMPOSITE_MULTIPLICATIVE) {
624: /* 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 */
625: if (!jac->Afield) {
626: PetscCalloc1(nsplit,&jac->Afield);
627: if (jac->offdiag_use_amat) {
628: MatGetSubMatrix(pc->mat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);
629: } else {
630: MatGetSubMatrix(pc->pmat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);
631: }
632: } else {
633: if (jac->offdiag_use_amat) {
634: MatGetSubMatrix(pc->mat,ilink->next->is,ilink->is,MAT_REUSE_MATRIX,&jac->Afield[1]);
635: } else {
636: MatGetSubMatrix(pc->pmat,ilink->next->is,ilink->is,MAT_REUSE_MATRIX,&jac->Afield[1]);
637: }
638: }
639: } else {
640: if (!jac->Afield) {
641: PetscMalloc1(nsplit,&jac->Afield);
642: for (i=0; i<nsplit; i++) {
643: if (jac->offdiag_use_amat) {
644: MatGetSubMatrix(pc->mat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);
645: } else {
646: MatGetSubMatrix(pc->pmat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);
647: }
648: ilink = ilink->next;
649: }
650: } else {
651: for (i=0; i<nsplit; i++) {
652: if (jac->offdiag_use_amat) {
653: MatGetSubMatrix(pc->mat,ilink->is,NULL,MAT_REUSE_MATRIX,&jac->Afield[i]);
654: } else {
655: MatGetSubMatrix(pc->pmat,ilink->is,NULL,MAT_REUSE_MATRIX,&jac->Afield[i]);
656: }
657: ilink = ilink->next;
658: }
659: }
660: }
661: }
663: if (jac->type == PC_COMPOSITE_SCHUR) {
664: IS ccis;
665: PetscInt rstart,rend;
666: char lscname[256];
667: PetscObject LSC_L;
669: if (nsplit != 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"To use Schur complement preconditioner you must have exactly 2 fields");
671: /* When extracting off-diagonal submatrices, we take complements from this range */
672: MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);
674: /* need to handle case when one is resetting up the preconditioner */
675: if (jac->schur) {
676: KSP kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper;
678: MatSchurComplementGetKSP(jac->schur, &kspInner);
679: ilink = jac->head;
680: ISComplement(ilink->is_col,rstart,rend,&ccis);
681: if (jac->offdiag_use_amat) {
682: MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);
683: } else {
684: MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);
685: }
686: ISDestroy(&ccis);
687: ilink = ilink->next;
688: ISComplement(ilink->is_col,rstart,rend,&ccis);
689: if (jac->offdiag_use_amat) {
690: MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);
691: } else {
692: MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);
693: }
694: ISDestroy(&ccis);
695: MatSchurComplementUpdateSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);
696: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
697: MatDestroy(&jac->schurp);
698: MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);
699: }
700: if (kspA != kspInner) {
701: KSPSetOperators(kspA,jac->mat[0],jac->pmat[0]);
702: }
703: if (kspUpper != kspA) {
704: KSPSetOperators(kspUpper,jac->mat[0],jac->pmat[0]);
705: }
706: KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));
707: } else {
708: const char *Dprefix;
709: char schurprefix[256], schurmatprefix[256];
710: char schurtestoption[256];
711: MatNullSpace sp;
712: PetscBool flg;
714: /* extract the A01 and A10 matrices */
715: ilink = jac->head;
716: ISComplement(ilink->is_col,rstart,rend,&ccis);
717: if (jac->offdiag_use_amat) {
718: MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
719: } else {
720: MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
721: }
722: ISDestroy(&ccis);
723: ilink = ilink->next;
724: ISComplement(ilink->is_col,rstart,rend,&ccis);
725: if (jac->offdiag_use_amat) {
726: MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
727: } else {
728: MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
729: }
730: ISDestroy(&ccis);
732: /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */
733: MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur);
734: MatSetType(jac->schur,MATSCHURCOMPLEMENT);
735: MatSchurComplementSetSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);
736: PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
737: /* 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? */
738: MatSetOptionsPrefix(jac->schur,schurmatprefix);
739: MatSetFromOptions(jac->schur);
740: MatGetNullSpace(jac->mat[1], &sp);
741: if (sp) {
742: MatSetNullSpace(jac->schur, sp);
743: }
745: PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);
746: PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
747: if (flg) {
748: DM dmInner;
749: KSP kspInner;
751: MatSchurComplementGetKSP(jac->schur, &kspInner);
752: PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
753: /* Indent this deeper to emphasize the "inner" nature of this solver. */
754: PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject) pc, 2);
755: KSPSetOptionsPrefix(kspInner, schurprefix);
757: /* Set DM for new solver */
758: KSPGetDM(jac->head->ksp, &dmInner);
759: KSPSetDM(kspInner, dmInner);
760: KSPSetDMActive(kspInner, PETSC_FALSE);
761: } else {
762: /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or
763: * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact,
764: * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for
765: * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make
766: * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used
767: * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */
768: KSPSetType(jac->head->ksp,KSPGMRES);
769: MatSchurComplementSetKSP(jac->schur,jac->head->ksp);
770: }
771: KSPSetOperators(jac->head->ksp,jac->mat[0],jac->pmat[0]);
772: KSPSetFromOptions(jac->head->ksp);
773: MatSetFromOptions(jac->schur);
775: PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);
776: PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
777: if (flg) {
778: DM dmInner;
780: PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
781: KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper);
782: KSPSetErrorIfNotConverged(jac->kspupper,pc->erroriffailure);
783: KSPSetOptionsPrefix(jac->kspupper, schurprefix);
784: KSPGetDM(jac->head->ksp, &dmInner);
785: KSPSetDM(jac->kspupper, dmInner);
786: KSPSetDMActive(jac->kspupper, PETSC_FALSE);
787: KSPSetFromOptions(jac->kspupper);
788: KSPSetOperators(jac->kspupper,jac->mat[0],jac->pmat[0]);
789: VecDuplicate(jac->head->x, &jac->head->z);
790: } else {
791: jac->kspupper = jac->head->ksp;
792: PetscObjectReference((PetscObject) jac->head->ksp);
793: }
795: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
796: MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);
797: }
798: KSPCreate(PetscObjectComm((PetscObject)pc),&jac->kspschur);
799: KSPSetErrorIfNotConverged(jac->kspschur,pc->erroriffailure);
800: PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);
801: PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);
802: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) {
803: PC pcschur;
804: KSPGetPC(jac->kspschur,&pcschur);
805: PCSetType(pcschur,PCNONE);
806: /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */
807: } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) {
808: MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user);
809: }
810: KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));
811: KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix);
812: KSPSetOptionsPrefix(jac->kspschur, Dprefix);
813: /* propogate DM */
814: {
815: DM sdm;
816: KSPGetDM(jac->head->next->ksp, &sdm);
817: if (sdm) {
818: KSPSetDM(jac->kspschur, sdm);
819: KSPSetDMActive(jac->kspschur, PETSC_FALSE);
820: }
821: }
822: /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
823: /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */
824: KSPSetFromOptions(jac->kspschur);
825: }
827: /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */
828: PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_L",ilink->splitname);
829: PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);
830: if (!LSC_L) {PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);}
831: if (LSC_L) {PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L);}
832: PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_Lp",ilink->splitname);
833: PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);
834: if (!LSC_L) {PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);}
835: if (LSC_L) {PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L);}
836: } else {
837: /* set up the individual splits' PCs */
838: i = 0;
839: ilink = jac->head;
840: while (ilink) {
841: KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i]);
842: /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
843: if (!jac->suboptionsset) {KSPSetFromOptions(ilink->ksp);}
844: i++;
845: ilink = ilink->next;
846: }
847: }
849: jac->suboptionsset = PETSC_TRUE;
850: return(0);
851: }
853: #define FieldSplitSplitSolveAdd(ilink,xx,yy) \
854: (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
855: VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
856: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\
857: KSPSolve(ilink->ksp,ilink->x,ilink->y) || \
858: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\
859: VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \
860: VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE))
864: static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y)
865: {
866: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
867: PetscErrorCode ierr;
868: PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
869: KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;
872: switch (jac->schurfactorization) {
873: case PC_FIELDSPLIT_SCHUR_FACT_DIAG:
874: /* [A00 0; 0 -S], positive definite, suitable for MINRES */
875: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
876: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
877: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
878: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
879: KSPSolve(kspA,ilinkA->x,ilinkA->y);
880: PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
881: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
882: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
883: PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
884: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
885: PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
886: VecScale(ilinkD->y,-1.);
887: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
888: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
889: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
890: break;
891: case PC_FIELDSPLIT_SCHUR_FACT_LOWER:
892: /* [A00 0; A10 S], suitable for left preconditioning */
893: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
894: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
895: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
896: KSPSolve(kspA,ilinkA->x,ilinkA->y);
897: PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
898: MatMult(jac->C,ilinkA->y,ilinkD->x);
899: VecScale(ilinkD->x,-1.);
900: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
901: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
902: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
903: PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
904: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
905: PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
906: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
907: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
908: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
909: break;
910: case PC_FIELDSPLIT_SCHUR_FACT_UPPER:
911: /* [A00 A01; 0 S], suitable for right preconditioning */
912: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
913: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
914: PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
915: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
916: PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL); MatMult(jac->B,ilinkD->y,ilinkA->x);
917: VecScale(ilinkA->x,-1.);
918: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);
919: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
920: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);
921: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
922: KSPSolve(kspA,ilinkA->x,ilinkA->y);
923: PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
924: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
925: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
926: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
927: break;
928: case PC_FIELDSPLIT_SCHUR_FACT_FULL:
929: /* [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 */
930: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
931: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
932: PetscLogEventBegin(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);
933: KSPSolve(kspLower,ilinkA->x,ilinkA->y);
934: PetscLogEventEnd(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);
935: MatMult(jac->C,ilinkA->y,ilinkD->x);
936: VecScale(ilinkD->x,-1.0);
937: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
938: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
940: PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
941: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
942: PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
943: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
944: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
946: if (kspUpper == kspA) {
947: MatMult(jac->B,ilinkD->y,ilinkA->y);
948: VecAXPY(ilinkA->x,-1.0,ilinkA->y);
949: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
950: KSPSolve(kspA,ilinkA->x,ilinkA->y);
951: PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
952: } else {
953: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
954: KSPSolve(kspA,ilinkA->x,ilinkA->y);
955: PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
956: MatMult(jac->B,ilinkD->y,ilinkA->x);
957: PetscLogEventBegin(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);
958: KSPSolve(kspUpper,ilinkA->x,ilinkA->z);
959: PetscLogEventEnd(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);
960: VecAXPY(ilinkA->y,-1.0,ilinkA->z);
961: }
962: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
963: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
964: }
965: return(0);
966: }
970: static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y)
971: {
972: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
973: PetscErrorCode ierr;
974: PC_FieldSplitLink ilink = jac->head;
975: PetscInt cnt,bs;
976: KSPConvergedReason reason;
979: if (jac->type == PC_COMPOSITE_ADDITIVE) {
980: if (jac->defaultsplit) {
981: VecGetBlockSize(x,&bs);
982: 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);
983: VecGetBlockSize(y,&bs);
984: 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);
985: VecStrideGatherAll(x,jac->x,INSERT_VALUES);
986: while (ilink) {
987: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
988: KSPSolve(ilink->ksp,ilink->x,ilink->y);
989: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
990: KSPGetConvergedReason(ilink->ksp,&reason);
991: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
992: pc->failedreason = PC_SUBPC_ERROR;
993: }
994: ilink = ilink->next;
995: }
996: VecStrideScatterAll(jac->y,y,INSERT_VALUES);
997: } else {
998: VecSet(y,0.0);
999: while (ilink) {
1000: FieldSplitSplitSolveAdd(ilink,x,y);
1001: KSPGetConvergedReason(ilink->ksp,&reason);
1002: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1003: pc->failedreason = PC_SUBPC_ERROR;
1004: }
1005: ilink = ilink->next;
1006: }
1007: }
1008: } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) {
1009: VecSet(y,0.0);
1010: /* solve on first block for first block variables */
1011: VecScatterBegin(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);
1012: VecScatterEnd(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);
1013: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1014: KSPSolve(ilink->ksp,ilink->x,ilink->y);
1015: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1016: KSPGetConvergedReason(ilink->ksp,&reason);
1017: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1018: pc->failedreason = PC_SUBPC_ERROR;
1019: }
1020: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1021: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1023: /* compute the residual only onto second block variables using first block variables */
1024: MatMult(jac->Afield[1],ilink->y,ilink->next->x);
1025: ilink = ilink->next;
1026: VecScale(ilink->x,-1.0);
1027: VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1028: VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1030: /* solve on second block variables */
1031: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1032: KSPSolve(ilink->ksp,ilink->x,ilink->y);
1033: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1034: KSPGetConvergedReason(ilink->ksp,&reason);
1035: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1036: pc->failedreason = PC_SUBPC_ERROR;
1037: }
1038: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1039: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1040: } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1041: if (!jac->w1) {
1042: VecDuplicate(x,&jac->w1);
1043: VecDuplicate(x,&jac->w2);
1044: }
1045: VecSet(y,0.0);
1046: FieldSplitSplitSolveAdd(ilink,x,y);
1047: KSPGetConvergedReason(ilink->ksp,&reason);
1048: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1049: pc->failedreason = PC_SUBPC_ERROR;
1050: }
1051: cnt = 1;
1052: while (ilink->next) {
1053: ilink = ilink->next;
1054: /* compute the residual only over the part of the vector needed */
1055: MatMult(jac->Afield[cnt++],y,ilink->x);
1056: VecScale(ilink->x,-1.0);
1057: VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1058: VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1059: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1060: KSPSolve(ilink->ksp,ilink->x,ilink->y);
1061: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1062: KSPGetConvergedReason(ilink->ksp,&reason);
1063: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1064: pc->failedreason = PC_SUBPC_ERROR;
1065: }
1066: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1067: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1068: }
1069: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1070: cnt -= 2;
1071: while (ilink->previous) {
1072: ilink = ilink->previous;
1073: /* compute the residual only over the part of the vector needed */
1074: MatMult(jac->Afield[cnt--],y,ilink->x);
1075: VecScale(ilink->x,-1.0);
1076: VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1077: VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1078: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1079: KSPSolve(ilink->ksp,ilink->x,ilink->y);
1080: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1081: KSPGetConvergedReason(ilink->ksp,&reason);
1082: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1083: pc->failedreason = PC_SUBPC_ERROR;
1084: }
1085: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1086: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1087: }
1088: }
1089: } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type);
1090: return(0);
1091: }
1093: #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \
1094: (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
1095: VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
1096: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \
1097: KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \
1098: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \
1099: VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \
1100: VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE))
1104: static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y)
1105: {
1106: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1107: PetscErrorCode ierr;
1108: PC_FieldSplitLink ilink = jac->head;
1109: PetscInt bs;
1110: KSPConvergedReason reason;
1113: if (jac->type == PC_COMPOSITE_ADDITIVE) {
1114: if (jac->defaultsplit) {
1115: VecGetBlockSize(x,&bs);
1116: 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);
1117: VecGetBlockSize(y,&bs);
1118: 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);
1119: VecStrideGatherAll(x,jac->x,INSERT_VALUES);
1120: while (ilink) {
1121: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1122: KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);
1123: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1124: KSPGetConvergedReason(ilink->ksp,&reason);
1125: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1126: pc->failedreason = PC_SUBPC_ERROR;
1127: }
1128: ilink = ilink->next;
1129: }
1130: VecStrideScatterAll(jac->y,y,INSERT_VALUES);
1131: } else {
1132: VecSet(y,0.0);
1133: while (ilink) {
1134: FieldSplitSplitSolveAddTranspose(ilink,x,y);
1135: KSPGetConvergedReason(ilink->ksp,&reason);
1136: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1137: pc->failedreason = PC_SUBPC_ERROR;
1138: }
1139: ilink = ilink->next;
1140: }
1141: }
1142: } else {
1143: if (!jac->w1) {
1144: VecDuplicate(x,&jac->w1);
1145: VecDuplicate(x,&jac->w2);
1146: }
1147: VecSet(y,0.0);
1148: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1149: FieldSplitSplitSolveAddTranspose(ilink,x,y);
1150: KSPGetConvergedReason(ilink->ksp,&reason);
1151: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1152: pc->failedreason = PC_SUBPC_ERROR;
1153: }
1154: while (ilink->next) {
1155: ilink = ilink->next;
1156: MatMultTranspose(pc->mat,y,jac->w1);
1157: VecWAXPY(jac->w2,-1.0,jac->w1,x);
1158: FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1159: }
1160: while (ilink->previous) {
1161: ilink = ilink->previous;
1162: MatMultTranspose(pc->mat,y,jac->w1);
1163: VecWAXPY(jac->w2,-1.0,jac->w1,x);
1164: FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1165: }
1166: } else {
1167: while (ilink->next) { /* get to last entry in linked list */
1168: ilink = ilink->next;
1169: }
1170: FieldSplitSplitSolveAddTranspose(ilink,x,y);
1171: KSPGetConvergedReason(ilink->ksp,&reason);
1172: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1173: pc->failedreason = PC_SUBPC_ERROR;
1174: }
1175: while (ilink->previous) {
1176: ilink = ilink->previous;
1177: MatMultTranspose(pc->mat,y,jac->w1);
1178: VecWAXPY(jac->w2,-1.0,jac->w1,x);
1179: FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1180: }
1181: }
1182: }
1183: return(0);
1184: }
1188: static PetscErrorCode PCReset_FieldSplit(PC pc)
1189: {
1190: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1191: PetscErrorCode ierr;
1192: PC_FieldSplitLink ilink = jac->head,next;
1195: while (ilink) {
1196: KSPReset(ilink->ksp);
1197: VecDestroy(&ilink->x);
1198: VecDestroy(&ilink->y);
1199: VecDestroy(&ilink->z);
1200: VecScatterDestroy(&ilink->sctx);
1201: if (!ilink->is_orig) { /* save the original IS */
1202: PetscObjectReference((PetscObject)ilink->is);
1203: ilink->is_orig = ilink->is;
1204: }
1205: ISDestroy(&ilink->is);
1206: ISDestroy(&ilink->is_col);
1207: next = ilink->next;
1208: ilink = next;
1209: }
1210: PetscFree2(jac->x,jac->y);
1211: if (jac->mat && jac->mat != jac->pmat) {
1212: MatDestroyMatrices(jac->nsplits,&jac->mat);
1213: } else if (jac->mat) {
1214: jac->mat = NULL;
1215: }
1216: if (jac->pmat) {MatDestroyMatrices(jac->nsplits,&jac->pmat);}
1217: if (jac->Afield) {MatDestroyMatrices(jac->nsplits,&jac->Afield);}
1218: VecDestroy(&jac->w1);
1219: VecDestroy(&jac->w2);
1220: MatDestroy(&jac->schur);
1221: MatDestroy(&jac->schurp);
1222: MatDestroy(&jac->schur_user);
1223: KSPDestroy(&jac->kspschur);
1224: KSPDestroy(&jac->kspupper);
1225: MatDestroy(&jac->B);
1226: MatDestroy(&jac->C);
1227: jac->reset = PETSC_TRUE;
1228: jac->isrestrict = PETSC_FALSE;
1229: return(0);
1230: }
1234: static PetscErrorCode PCDestroy_FieldSplit(PC pc)
1235: {
1236: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1237: PetscErrorCode ierr;
1238: PC_FieldSplitLink ilink = jac->head,next;
1241: PCReset_FieldSplit(pc);
1242: while (ilink) {
1243: KSPDestroy(&ilink->ksp);
1244: ISDestroy(&ilink->is_orig);
1245: next = ilink->next;
1246: PetscFree(ilink->splitname);
1247: PetscFree(ilink->fields);
1248: PetscFree(ilink->fields_col);
1249: PetscFree(ilink);
1250: ilink = next;
1251: }
1252: PetscFree2(jac->x,jac->y);
1253: PetscFree(pc->data);
1254: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL);
1255: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL);
1256: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL);
1257: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL);
1258: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL);
1259: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL);
1260: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL);
1261: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL);
1262: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",NULL);
1263: return(0);
1264: }
1268: static PetscErrorCode PCSetFromOptions_FieldSplit(PetscOptionItems *PetscOptionsObject,PC pc)
1269: {
1270: PetscErrorCode ierr;
1271: PetscInt bs;
1272: PetscBool flg,stokes = PETSC_FALSE;
1273: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1274: PCCompositeType ctype;
1277: PetscOptionsHead(PetscOptionsObject,"FieldSplit options");
1278: PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL);
1279: PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);
1280: if (flg) {
1281: PCFieldSplitSetBlockSize(pc,bs);
1282: }
1283: jac->diag_use_amat = pc->useAmat;
1284: 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);
1285: jac->offdiag_use_amat = pc->useAmat;
1286: 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);
1287: /* FIXME: No programmatic equivalent to the following. */
1288: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);
1289: if (stokes) {
1290: PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);
1291: jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF;
1292: }
1294: PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);
1295: if (flg) {
1296: PCFieldSplitSetType(pc,ctype);
1297: }
1298: /* Only setup fields once */
1299: if ((jac->bs > 0) && (jac->nsplits == 0)) {
1300: /* only allow user to set fields from command line if bs is already known.
1301: otherwise user can set them in PCFieldSplitSetDefaults() */
1302: PCFieldSplitSetRuntimeSplits_Private(pc);
1303: if (jac->splitdefined) {PetscInfo(pc,"Splits defined using the options database\n");}
1304: }
1305: if (jac->type == PC_COMPOSITE_SCHUR) {
1306: PetscOptionsGetEnum(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);
1307: if (flg) {PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");}
1308: 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);
1309: PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSetSchurPre",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL);
1310: }
1311: PetscOptionsTail();
1312: return(0);
1313: }
1315: /*------------------------------------------------------------------------------------*/
1319: static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1320: {
1321: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1322: PetscErrorCode ierr;
1323: PC_FieldSplitLink ilink,next = jac->head;
1324: char prefix[128];
1325: PetscInt i;
1328: if (jac->splitdefined) {
1329: PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1330: return(0);
1331: }
1332: for (i=0; i<n; i++) {
1333: if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs);
1334: if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]);
1335: }
1336: PetscNew(&ilink);
1337: if (splitname) {
1338: PetscStrallocpy(splitname,&ilink->splitname);
1339: } else {
1340: PetscMalloc1(3,&ilink->splitname);
1341: PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);
1342: }
1343: 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 */
1344: PetscMalloc1(n,&ilink->fields);
1345: PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));
1346: PetscMalloc1(n,&ilink->fields_col);
1347: PetscMemcpy(ilink->fields_col,fields_col,n*sizeof(PetscInt));
1349: ilink->nfields = n;
1350: ilink->next = NULL;
1351: KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1352: KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);
1353: PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1354: KSPSetType(ilink->ksp,KSPPREONLY);
1355: PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);
1357: PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);
1358: KSPSetOptionsPrefix(ilink->ksp,prefix);
1360: if (!next) {
1361: jac->head = ilink;
1362: ilink->previous = NULL;
1363: } else {
1364: while (next->next) {
1365: next = next->next;
1366: }
1367: next->next = ilink;
1368: ilink->previous = next;
1369: }
1370: jac->nsplits++;
1371: return(0);
1372: }
1376: static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp)
1377: {
1378: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1382: PetscMalloc1(jac->nsplits,subksp);
1383: MatSchurComplementGetKSP(jac->schur,*subksp);
1385: (*subksp)[1] = jac->kspschur;
1386: if (n) *n = jac->nsplits;
1387: return(0);
1388: }
1392: static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp)
1393: {
1394: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1395: PetscErrorCode ierr;
1396: PetscInt cnt = 0;
1397: PC_FieldSplitLink ilink = jac->head;
1400: PetscMalloc1(jac->nsplits,subksp);
1401: while (ilink) {
1402: (*subksp)[cnt++] = ilink->ksp;
1403: ilink = ilink->next;
1404: }
1405: 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);
1406: if (n) *n = jac->nsplits;
1407: return(0);
1408: }
1412: /*@C
1413: PCFieldSplitRestrictIS - Restricts the fieldsplit ISs to be within a given IS.
1415: Input Parameters:
1416: + pc - the preconditioner context
1417: + is - the index set that defines the indices to which the fieldsplit is to be restricted
1419: Level: advanced
1421: @*/
1422: PetscErrorCode PCFieldSplitRestrictIS(PC pc,IS isy)
1423: {
1429: PetscTryMethod(pc,"PCFieldSplitRestrictIS_C",(PC,IS),(pc,isy));
1430: return(0);
1431: }
1436: static PetscErrorCode PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy)
1437: {
1438: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1439: PetscErrorCode ierr;
1440: PC_FieldSplitLink ilink = jac->head, next;
1441: PetscInt localsize,size,sizez,i;
1442: const PetscInt *ind, *indz;
1443: PetscInt *indc, *indcz;
1444: PetscBool flg;
1447: ISGetLocalSize(isy,&localsize);
1448: MPI_Scan(&localsize,&size,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)isy));
1449: size -= localsize;
1450: while(ilink) {
1451: IS isrl,isr;
1452: PC subpc;
1453: if (jac->reset) {
1454: ISEmbed(ilink->is_orig, isy, PETSC_TRUE, &isrl);
1455: } else {
1456: ISEmbed(ilink->is, isy, PETSC_TRUE, &isrl);
1457: }
1458: ISGetLocalSize(isrl,&localsize);
1459: PetscMalloc1(localsize,&indc);
1460: ISGetIndices(isrl,&ind);
1461: PetscMemcpy(indc,ind,localsize*sizeof(PetscInt));
1462: ISRestoreIndices(isrl,&ind);
1463: ISDestroy(&isrl);
1464: for (i=0; i<localsize; i++) *(indc+i) += size;
1465: ISCreateGeneral(PetscObjectComm((PetscObject)isy),localsize,indc,PETSC_OWN_POINTER,&isr);
1466: PetscObjectReference((PetscObject)isr);
1467: ISDestroy(&ilink->is);
1468: ilink->is = isr;
1469: PetscObjectReference((PetscObject)isr);
1470: ISDestroy(&ilink->is_col);
1471: ilink->is_col = isr;
1472: ISDestroy(&isr);
1473: KSPGetPC(ilink->ksp, &subpc);
1474: PetscObjectTypeCompare((PetscObject)subpc,PCFIELDSPLIT,&flg);
1475: if(flg) {
1476: IS iszl,isz;
1477: MPI_Comm comm;
1478: if (jac->reset) {
1479: ISGetLocalSize(ilink->is_orig,&localsize);
1480: comm = PetscObjectComm((PetscObject)ilink->is_orig);
1481: ISEmbed(isy, ilink->is_orig, PETSC_TRUE, &iszl);
1482: } else {
1483: ISGetLocalSize(ilink->is,&localsize);
1484: comm = PetscObjectComm((PetscObject)ilink->is);
1485: ISEmbed(isy, ilink->is, PETSC_TRUE, &iszl);
1486: }
1487: MPI_Scan(&localsize,&sizez,1,MPIU_INT,MPI_SUM,comm);
1488: sizez -= localsize;
1489: ISGetLocalSize(iszl,&localsize);
1490: PetscMalloc1(localsize,&indcz);
1491: ISGetIndices(iszl,&indz);
1492: PetscMemcpy(indcz,indz,localsize*sizeof(PetscInt));
1493: ISRestoreIndices(iszl,&indz);
1494: ISDestroy(&iszl);
1495: for (i=0; i<localsize; i++) *(indcz+i) += sizez;
1496: ISCreateGeneral(comm,localsize,indcz,PETSC_OWN_POINTER,&isz);
1497: PCFieldSplitRestrictIS(subpc,isz);
1498: ISDestroy(&isz);
1499: }
1500: next = ilink->next;
1501: ilink = next;
1502: }
1503: jac->isrestrict = PETSC_TRUE;
1504: return(0);
1505: }
1509: static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is)
1510: {
1511: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1512: PetscErrorCode ierr;
1513: PC_FieldSplitLink ilink, next = jac->head;
1514: char prefix[128];
1517: if (jac->splitdefined) {
1518: PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1519: return(0);
1520: }
1521: PetscNew(&ilink);
1522: if (splitname) {
1523: PetscStrallocpy(splitname,&ilink->splitname);
1524: } else {
1525: PetscMalloc1(8,&ilink->splitname);
1526: PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);
1527: }
1528: 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 */
1529: PetscObjectReference((PetscObject)is);
1530: ISDestroy(&ilink->is);
1531: ilink->is = is;
1532: PetscObjectReference((PetscObject)is);
1533: ISDestroy(&ilink->is_col);
1534: ilink->is_col = is;
1535: ilink->next = NULL;
1536: KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1537: KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);
1538: PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1539: KSPSetType(ilink->ksp,KSPPREONLY);
1540: PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);
1542: PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);
1543: KSPSetOptionsPrefix(ilink->ksp,prefix);
1545: if (!next) {
1546: jac->head = ilink;
1547: ilink->previous = NULL;
1548: } else {
1549: while (next->next) {
1550: next = next->next;
1551: }
1552: next->next = ilink;
1553: ilink->previous = next;
1554: }
1555: jac->nsplits++;
1556: return(0);
1557: }
1561: /*@
1562: PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner
1564: Logically Collective on PC
1566: Input Parameters:
1567: + pc - the preconditioner context
1568: . splitname - name of this split, if NULL the number of the split is used
1569: . n - the number of fields in this split
1570: - fields - the fields in this split
1572: Level: intermediate
1574: Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field.
1576: The PCFieldSplitSetFields() is for defining fields as strided blocks. For example, if the block
1577: 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
1578: 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x....
1579: where the numbered entries indicate what is in the field.
1581: This function is called once per split (it creates a new split each time). Solve options
1582: for this split will be available under the prefix -fieldsplit_SPLITNAME_.
1584: Developer Note: This routine does not actually create the IS representing the split, that is delayed
1585: until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be
1586: available when this routine is called.
1588: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS()
1590: @*/
1591: PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1592: {
1598: if (n < 1) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname);
1600: PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt*,const PetscInt*),(pc,splitname,n,fields,fields_col));
1601: return(0);
1602: }
1606: /*@
1607: PCFieldSplitSetDiagUseAmat - set flag indicating whether to extract diagonal blocks from Amat (rather than Pmat)
1609: Logically Collective on PC
1611: Input Parameters:
1612: + pc - the preconditioner object
1613: - flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from
1615: Options Database:
1616: . -pc_fieldsplit_diag_use_amat
1618: Level: intermediate
1620: .seealso: PCFieldSplitGetDiagUseAmat(), PCFieldSplitSetOffDiagUseAmat(), PCFIELDSPLIT
1622: @*/
1623: PetscErrorCode PCFieldSplitSetDiagUseAmat(PC pc,PetscBool flg)
1624: {
1625: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1626: PetscBool isfs;
1631: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1632: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1633: jac->diag_use_amat = flg;
1634: return(0);
1635: }
1639: /*@
1640: PCFieldSplitGetDiagUseAmat - get the flag indicating whether to extract 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 diagonal blocks from
1651: Level: intermediate
1653: .seealso: PCFieldSplitSetDiagUseAmat(), PCFieldSplitGetOffDiagUseAmat(), PCFIELDSPLIT
1655: @*/
1656: PetscErrorCode PCFieldSplitGetDiagUseAmat(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->diag_use_amat;
1668: return(0);
1669: }
1673: /*@
1674: PCFieldSplitSetOffDiagUseAmat - set flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat)
1676: Logically Collective on PC
1678: Input Parameters:
1679: + pc - the preconditioner object
1680: - flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from
1682: Options Database:
1683: . -pc_fieldsplit_off_diag_use_amat
1685: Level: intermediate
1687: .seealso: PCFieldSplitGetOffDiagUseAmat(), PCFieldSplitSetDiagUseAmat(), PCFIELDSPLIT
1689: @*/
1690: PetscErrorCode PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg)
1691: {
1692: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1693: PetscBool isfs;
1698: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1699: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1700: jac->offdiag_use_amat = flg;
1701: return(0);
1702: }
1706: /*@
1707: PCFieldSplitGetOffDiagUseAmat - get the flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat)
1709: Logically Collective on PC
1711: Input Parameters:
1712: . pc - the preconditioner object
1714: Output Parameters:
1715: . flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from
1718: Level: intermediate
1720: .seealso: PCFieldSplitSetOffDiagUseAmat(), PCFieldSplitGetDiagUseAmat(), PCFIELDSPLIT
1722: @*/
1723: PetscErrorCode PCFieldSplitGetOffDiagUseAmat(PC pc,PetscBool *flg)
1724: {
1725: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1726: PetscBool isfs;
1732: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1733: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1734: *flg = jac->offdiag_use_amat;
1735: return(0);
1736: }
1742: /*@C
1743: PCFieldSplitSetIS - Sets the exact elements for field
1745: Logically Collective on PC
1747: Input Parameters:
1748: + pc - the preconditioner context
1749: . splitname - name of this split, if NULL the number of the split is used
1750: - is - the index set that defines the vector elements in this field
1753: Notes:
1754: Use PCFieldSplitSetFields(), for fields defined by strided types.
1756: This function is called once per split (it creates a new split each time). Solve options
1757: for this split will be available under the prefix -fieldsplit_SPLITNAME_.
1759: Level: intermediate
1761: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize()
1763: @*/
1764: PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is)
1765: {
1772: PetscUseMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));
1773: return(0);
1774: }
1778: /*@
1779: PCFieldSplitGetIS - Retrieves the elements for a field as an IS
1781: Logically Collective on PC
1783: Input Parameters:
1784: + pc - the preconditioner context
1785: - splitname - name of this split
1787: Output Parameter:
1788: - is - the index set that defines the vector elements in this field, or NULL if the field is not found
1790: Level: intermediate
1792: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS()
1794: @*/
1795: PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is)
1796: {
1803: {
1804: PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;
1805: PC_FieldSplitLink ilink = jac->head;
1806: PetscBool found;
1808: *is = NULL;
1809: while (ilink) {
1810: PetscStrcmp(ilink->splitname, splitname, &found);
1811: if (found) {
1812: *is = ilink->is;
1813: break;
1814: }
1815: ilink = ilink->next;
1816: }
1817: }
1818: return(0);
1819: }
1823: /*@
1824: PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the
1825: fieldsplit preconditioner. If not set the matrix block size is used.
1827: Logically Collective on PC
1829: Input Parameters:
1830: + pc - the preconditioner context
1831: - bs - the block size
1833: Level: intermediate
1835: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields()
1837: @*/
1838: PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs)
1839: {
1845: PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));
1846: return(0);
1847: }
1851: /*@C
1852: PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits
1854: Collective on KSP
1856: Input Parameter:
1857: . pc - the preconditioner context
1859: Output Parameters:
1860: + n - the number of splits
1861: - pc - the array of KSP contexts
1863: Note:
1864: After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user
1865: (not the KSP just the array that contains them).
1867: You must call KSPSetUp() before calling PCFieldSplitGetSubKSP().
1869: Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs.
1870: You can call PCFieldSplitGetSubKSP(pc,n,NULL_OBJECT,ierr) to determine how large the
1871: KSP array must be.
1874: Level: advanced
1876: .seealso: PCFIELDSPLIT
1877: @*/
1878: PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
1879: {
1885: PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));
1886: return(0);
1887: }
1891: /*@
1892: PCFieldSplitSetSchurPre - Indicates if the Schur complement is preconditioned by a preconditioner constructed by the
1893: A11 matrix. Otherwise no preconditioner is used.
1895: Collective on PC
1897: Input Parameters:
1898: + pc - the preconditioner context
1899: . 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
1900: PC_FIELDSPLIT_SCHUR_PRE_SELFP, and PC_FIELDSPLIT_SCHUR_PRE_FULL
1901: - userpre - matrix to use for preconditioning, or NULL
1903: Options Database:
1904: . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11. See notes for meaning of various arguments
1906: Notes:
1907: $ If ptype is
1908: $ a11 then the preconditioner for the Schur complement is generated from the block diagonal part of the preconditioner
1909: $ matrix associated with the Schur complement (i.e. A11), not he Schur complement matrix
1910: $ self the preconditioner for the Schur complement is generated from the symbolic representation of the Schur complement matrix:
1911: $ The only preconditioner that currently works with this symbolic respresentation matrix object is the PCLSC
1912: $ preconditioner
1913: $ user then the preconditioner for the Schur complement is generated from the user provided matrix (pre argument
1914: $ to this function).
1915: $ selfp then the preconditioning for the Schur complement is generated from an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01
1916: $ This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be
1917: $ lumped before extracting the diagonal using the additional option -fieldsplit_1_mat_schur_complement_ainv_type lump
1918: $ full then the preconditioner for the Schur complement is generated from the exact Schur complement matrix representation computed internally by PFIELDSPLIT (this is expensive)
1919: $ useful mostly as a test that the Schur complement approach can work for your problem
1921: When solving a saddle point problem, where the A11 block is identically zero, using a11 as the ptype only makes sense
1922: with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and
1923: -fieldsplit_1_pc_type lsc which uses the least squares commutator to compute a preconditioner for the Schur complement.
1925: Level: intermediate
1927: .seealso: PCFieldSplitGetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType,
1928: MatSchurComplementSetAinvType(), PCLSC
1930: @*/
1931: PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
1932: {
1937: PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));
1938: return(0);
1939: }
1940: PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) {return PCFieldSplitSetSchurPre(pc,ptype,pre);} /* Deprecated name */
1944: /*@
1945: PCFieldSplitGetSchurPre - For Schur complement fieldsplit, determine how the Schur complement will be
1946: preconditioned. See PCFieldSplitSetSchurPre() for details.
1948: Logically Collective on PC
1950: Input Parameters:
1951: . pc - the preconditioner context
1953: Output Parameters:
1954: + ptype - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11, PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_PRE_USER
1955: - userpre - matrix to use for preconditioning (with PC_FIELDSPLIT_PRE_USER), or NULL
1957: Level: intermediate
1959: .seealso: PCFieldSplitSetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC
1961: @*/
1962: PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
1963: {
1968: PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre));
1969: return(0);
1970: }
1974: /*@
1975: PCFieldSplitSchurGetS - extract the MatSchurComplement object used by this PC in case it needs to be configured separately
1977: Not collective
1979: Input Parameter:
1980: . pc - the preconditioner context
1982: Output Parameter:
1983: . S - the Schur complement matrix
1985: Notes:
1986: This matrix should not be destroyed using MatDestroy(); rather, use PCFieldSplitSchurRestoreS().
1988: Level: advanced
1990: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurRestoreS()
1992: @*/
1993: PetscErrorCode PCFieldSplitSchurGetS(PC pc,Mat *S)
1994: {
1996: const char* t;
1997: PetscBool isfs;
1998: PC_FieldSplit *jac;
2002: PetscObjectGetType((PetscObject)pc,&t);
2003: PetscStrcmp(t,PCFIELDSPLIT,&isfs);
2004: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
2005: jac = (PC_FieldSplit*)pc->data;
2006: if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
2007: if (S) *S = jac->schur;
2008: return(0);
2009: }
2013: /*@
2014: PCFieldSplitSchurRestoreS - restores the MatSchurComplement object used by this PC
2016: Not collective
2018: Input Parameters:
2019: + pc - the preconditioner context
2020: . S - the Schur complement matrix
2022: Level: advanced
2024: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurGetS()
2026: @*/
2027: PetscErrorCode PCFieldSplitSchurRestoreS(PC pc,Mat *S)
2028: {
2030: const char* t;
2031: PetscBool isfs;
2032: PC_FieldSplit *jac;
2036: PetscObjectGetType((PetscObject)pc,&t);
2037: PetscStrcmp(t,PCFIELDSPLIT,&isfs);
2038: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
2039: jac = (PC_FieldSplit*)pc->data;
2040: if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
2041: if (!S || *S != jac->schur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten");
2042: return(0);
2043: }
2048: static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
2049: {
2050: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2054: jac->schurpre = ptype;
2055: if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) {
2056: MatDestroy(&jac->schur_user);
2057: jac->schur_user = pre;
2058: PetscObjectReference((PetscObject)jac->schur_user);
2059: }
2060: return(0);
2061: }
2065: static PetscErrorCode PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
2066: {
2067: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2070: *ptype = jac->schurpre;
2071: *pre = jac->schur_user;
2072: return(0);
2073: }
2077: /*@
2078: PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain
2080: Collective on PC
2082: Input Parameters:
2083: + pc - the preconditioner context
2084: - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default
2086: Options Database:
2087: . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full
2090: Level: intermediate
2092: Notes:
2093: The FULL factorization is
2095: $ (A B) = (1 0) (A 0) (1 Ainv*B)
2096: $ (C D) (C*Ainv 1) (0 S) (0 1 )
2098: where S = D - 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,
2099: 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).
2101: If applied exactly, FULL factorization is a direct solver. The preconditioned operator with LOWER or UPPER has all eigenvalues equal to 1 and minimal polynomial
2102: of degree 2, so KSPGMRES converges in 2 iterations. If the iteration count is very low, consider using KSPFGMRES or KSPGCR which can use one less preconditioner
2103: application in this case. Note that the preconditioned operator may be highly non-normal, so such fast convergence may not be observed in practice. With DIAG,
2104: the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so KSPGMRES converges in at most 4 iterations.
2106: For symmetric problems in which A is positive definite and S is negative definite, DIAG can be used with KSPMINRES. Note that a flexible method like KSPFGMRES
2107: or KSPGCR must be used if the fieldsplit preconditioner is nonlinear (e.g. a few iterations of a Krylov method is used inside a split).
2109: References:
2110: + 1. - Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000).
2111: - 2. - Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001).
2113: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType
2114: @*/
2115: PetscErrorCode PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype)
2116: {
2121: PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));
2122: return(0);
2123: }
2127: static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype)
2128: {
2129: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2132: jac->schurfactorization = ftype;
2133: return(0);
2134: }
2138: /*@C
2139: PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement
2141: Collective on KSP
2143: Input Parameter:
2144: . pc - the preconditioner context
2146: Output Parameters:
2147: + A00 - the (0,0) block
2148: . A01 - the (0,1) block
2149: . A10 - the (1,0) block
2150: - A11 - the (1,1) block
2152: Level: advanced
2154: .seealso: PCFIELDSPLIT
2155: @*/
2156: PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11)
2157: {
2158: PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;
2162: if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach.");
2163: if (A00) *A00 = jac->pmat[0];
2164: if (A01) *A01 = jac->B;
2165: if (A10) *A10 = jac->C;
2166: if (A11) *A11 = jac->pmat[1];
2167: return(0);
2168: }
2172: static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type)
2173: {
2174: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2178: jac->type = type;
2179: if (type == PC_COMPOSITE_SCHUR) {
2180: pc->ops->apply = PCApply_FieldSplit_Schur;
2181: pc->ops->view = PCView_FieldSplit_Schur;
2183: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);
2184: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit);
2185: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit);
2186: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);
2188: } else {
2189: pc->ops->apply = PCApply_FieldSplit;
2190: pc->ops->view = PCView_FieldSplit;
2192: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2193: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",0);
2194: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",0);
2195: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);
2196: }
2197: return(0);
2198: }
2202: static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs)
2203: {
2204: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2207: if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs);
2208: 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);
2209: jac->bs = bs;
2210: return(0);
2211: }
2215: /*@
2216: PCFieldSplitSetType - Sets the type of fieldsplit preconditioner.
2218: Collective on PC
2220: Input Parameter:
2221: . pc - the preconditioner context
2222: . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR
2224: Options Database Key:
2225: . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type
2227: Level: Intermediate
2229: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
2231: .seealso: PCCompositeSetType()
2233: @*/
2234: PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type)
2235: {
2240: PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));
2241: return(0);
2242: }
2246: /*@
2247: PCFieldSplitGetType - Gets the type of fieldsplit preconditioner.
2249: Not collective
2251: Input Parameter:
2252: . pc - the preconditioner context
2254: Output Parameter:
2255: . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR
2257: Level: Intermediate
2259: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
2260: .seealso: PCCompositeSetType()
2261: @*/
2262: PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type)
2263: {
2264: PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;
2269: *type = jac->type;
2270: return(0);
2271: }
2275: /*@
2276: PCFieldSplitSetDMSplits - Flags whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible.
2278: Logically Collective
2280: Input Parameters:
2281: + pc - the preconditioner context
2282: - flg - boolean indicating whether to use field splits defined by the DM
2284: Options Database Key:
2285: . -pc_fieldsplit_dm_splits
2287: Level: Intermediate
2289: .keywords: PC, DM, composite preconditioner, additive, multiplicative
2291: .seealso: PCFieldSplitGetDMSplits()
2293: @*/
2294: PetscErrorCode PCFieldSplitSetDMSplits(PC pc,PetscBool flg)
2295: {
2296: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2297: PetscBool isfs;
2303: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2304: if (isfs) {
2305: jac->dm_splits = flg;
2306: }
2307: return(0);
2308: }
2313: /*@
2314: PCFieldSplitGetDMSplits - Returns flag indicating whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible.
2316: Logically Collective
2318: Input Parameter:
2319: . pc - the preconditioner context
2321: Output Parameter:
2322: . flg - boolean indicating whether to use field splits defined by the DM
2324: Level: Intermediate
2326: .keywords: PC, DM, composite preconditioner, additive, multiplicative
2328: .seealso: PCFieldSplitSetDMSplits()
2330: @*/
2331: PetscErrorCode PCFieldSplitGetDMSplits(PC pc,PetscBool* flg)
2332: {
2333: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2334: PetscBool isfs;
2340: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2341: if (isfs) {
2342: if(flg) *flg = jac->dm_splits;
2343: }
2344: return(0);
2345: }
2347: /* -------------------------------------------------------------------------------------*/
2348: /*MC
2349: PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual
2350: fields or groups of fields. See the users manual section "Solving Block Matrices" for more details.
2352: To set options on the solvers for each block append -fieldsplit_ to all the PC
2353: options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1
2355: To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP()
2356: and set the options directly on the resulting KSP object
2358: Level: intermediate
2360: Options Database Keys:
2361: + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split
2362: . -pc_fieldsplit_default - automatically add any fields to additional splits that have not
2363: been supplied explicitly by -pc_fieldsplit_%d_fields
2364: . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields)
2365: . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur> - type of relaxation or factorization splitting
2366: . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11; see PCFieldSplitSetSchurPre()
2367: . -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero or negative diagonal and uses Schur complement with no preconditioner as the solver
2369: - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_
2370: for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields
2372: Notes:
2373: Use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS()
2374: to define a field by an arbitrary collection of entries.
2376: If no fields are set the default is used. The fields are defined by entries strided by bs,
2377: beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(),
2378: if this is not called the block size defaults to the blocksize of the second matrix passed
2379: to KSPSetOperators()/PCSetOperators().
2381: $ For the Schur complement preconditioner if J = ( A00 A01 )
2382: $ ( A10 A11 )
2383: $ the preconditioner using full factorization is
2384: $ ( I -ksp(A00) A01 ) ( inv(A00) 0 ) ( I 0 )
2385: $ ( 0 I ) ( 0 ksp(S) ) ( -A10 ksp(A00) I )
2386: where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. S is the Schur complement
2387: $ S = A11 - A10 ksp(A00) A01
2388: 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
2389: in providing the SECOND split or 1 if not give). For PCFieldSplitGetKSP() when field number is 0,
2390: it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default
2391: A11 is used to construct a preconditioner for S, use PCFieldSplitSetSchurPre() for all the possible ways to construct the preconditioner for S.
2393: The factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above,
2394: diag gives
2395: $ ( inv(A00) 0 )
2396: $ ( 0 -ksp(S) )
2397: note that slightly counter intuitively there is a negative in front of the ksp(S) so that the preconditioner is positive definite. The lower factorization is the inverse of
2398: $ ( A00 0 )
2399: $ ( A10 S )
2400: where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of
2401: $ ( A00 A01 )
2402: $ ( 0 S )
2403: where again the inverses of A00 and S are applied using KSPs.
2405: If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS
2406: is used automatically for a second block.
2408: The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1.
2409: Generally it should be used with the AIJ format.
2411: The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see,
2412: for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT
2413: inside a smoother resulting in "Distributive Smoothers".
2415: Concepts: physics based preconditioners, block preconditioners
2417: There is a nice discussion of block preconditioners in
2419: [El08] A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier-Stokes equations
2420: Howard Elman, V.E. Howle, John Shadid, Robert Shuttleworth, Ray Tuminaro, Journal of Computational Physics 227 (2008) 1790--1808
2421: http://chess.cs.umd.edu/~elman/papers/tax.pdf
2423: The Constrained Pressure Preconditioner (CPR) does not appear to be currently implementable directly with PCFIELDSPLIT. CPR solves first the Schur complemented pressure equation, updates the
2424: residual on all variables and then applies a simple ILU like preconditioner on all the variables. So it is very much like the full Schur complement with selfp representing the Schur complement but instead
2425: of backsolving for the saturations in the last step it solves a full coupled (ILU) system for updates to all the variables.
2427: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC,
2428: PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre(),
2429: MatSchurComplementSetAinvType()
2430: M*/
2434: PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc)
2435: {
2437: PC_FieldSplit *jac;
2440: PetscNewLog(pc,&jac);
2442: jac->bs = -1;
2443: jac->nsplits = 0;
2444: jac->type = PC_COMPOSITE_MULTIPLICATIVE;
2445: jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */
2446: jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL;
2447: jac->dm_splits = PETSC_TRUE;
2449: pc->data = (void*)jac;
2451: pc->ops->apply = PCApply_FieldSplit;
2452: pc->ops->applytranspose = PCApplyTranspose_FieldSplit;
2453: pc->ops->setup = PCSetUp_FieldSplit;
2454: pc->ops->reset = PCReset_FieldSplit;
2455: pc->ops->destroy = PCDestroy_FieldSplit;
2456: pc->ops->setfromoptions = PCSetFromOptions_FieldSplit;
2457: pc->ops->view = PCView_FieldSplit;
2458: pc->ops->applyrichardson = 0;
2460: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2461: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);
2462: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);
2463: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);
2464: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);
2465: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",PCFieldSplitRestrictIS_FieldSplit);
2466: return(0);
2467: }