Actual source code: fieldsplit.c
petsc-3.10.5 2019-03-28
3: #include <petsc/private/pcimpl.h>
4: #include <petsc/private/kspimpl.h> /* This is needed to provide the appropriate PETSC_EXTERN for KSP_Solve_FS ....*/
5: #include <petscdm.h>
7: const char *const PCFieldSplitSchurPreTypes[] = {"SELF","SELFP","A11","USER","FULL","PCFieldSplitSchurPreType","PC_FIELDSPLIT_SCHUR_PRE_",0};
8: const char *const PCFieldSplitSchurFactTypes[] = {"DIAG","LOWER","UPPER","FULL","PCFieldSplitSchurFactType","PC_FIELDSPLIT_SCHUR_FACT_",0};
10: PetscLogEvent KSP_Solve_FS_0,KSP_Solve_FS_1,KSP_Solve_FS_S,KSP_Solve_FS_U,KSP_Solve_FS_L,KSP_Solve_FS_2,KSP_Solve_FS_3,KSP_Solve_FS_4;
12: typedef struct _PC_FieldSplitLink *PC_FieldSplitLink;
13: struct _PC_FieldSplitLink {
14: KSP ksp;
15: Vec x,y,z;
16: char *splitname;
17: PetscInt nfields;
18: PetscInt *fields,*fields_col;
19: VecScatter sctx;
20: IS is,is_col;
21: PC_FieldSplitLink next,previous;
22: PetscLogEvent event;
23: };
25: typedef struct {
26: PCCompositeType type;
27: PetscBool defaultsplit; /* Flag for a system with a set of 'k' scalar fields with the same layout (and bs = k) */
28: PetscBool splitdefined; /* Flag is set after the splits have been defined, to prevent more splits from being added */
29: PetscInt bs; /* Block size for IS and Mat structures */
30: PetscInt nsplits; /* Number of field divisions defined */
31: Vec *x,*y,w1,w2;
32: Mat *mat; /* The diagonal block for each split */
33: Mat *pmat; /* The preconditioning diagonal block for each split */
34: Mat *Afield; /* The rows of the matrix associated with each split */
35: PetscBool issetup;
37: /* Only used when Schur complement preconditioning is used */
38: Mat B; /* The (0,1) block */
39: Mat C; /* The (1,0) block */
40: Mat schur; /* The Schur complement S = A11 - A10 A00^{-1} A01, the KSP here, kspinner, is H_1 in [El08] */
41: Mat schurp; /* Assembled approximation to S built by MatSchurComplement to be used as a preconditioning matrix when solving with S */
42: Mat schur_user; /* User-provided preconditioning matrix for the Schur complement */
43: PCFieldSplitSchurPreType schurpre; /* Determines which preconditioning matrix is used for the Schur complement */
44: PCFieldSplitSchurFactType schurfactorization;
45: KSP kspschur; /* The solver for S */
46: KSP kspupper; /* The solver for A in the upper diagonal part of the factorization (H_2 in [El08]) */
47: PetscScalar schurscale; /* Scaling factor for the Schur complement solution with DIAG factorization */
49: PC_FieldSplitLink head;
50: PetscBool isrestrict; /* indicates PCFieldSplitRestrictIS() has been last called on this object, hack */
51: PetscBool suboptionsset; /* Indicates that the KSPSetFromOptions() has been called on the sub-KSPs */
52: PetscBool dm_splits; /* Whether to use DMCreateFieldDecomposition() whenever possible */
53: PetscBool diag_use_amat; /* Whether to extract diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */
54: PetscBool offdiag_use_amat; /* Whether to extract off-diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */
55: PetscBool detect; /* Whether to form 2-way split by finding zero diagonal entries */
56: } PC_FieldSplit;
58: /*
59: Notes:
60: there is no particular reason that pmat, x, and y are stored as arrays in PC_FieldSplit instead of
61: inside PC_FieldSplitLink, just historical. If you want to be able to add new fields after already using the
62: PC you could change this.
63: */
65: /* This helper is so that setting a user-provided preconditioning matrix is orthogonal to choosing to use it. This way the
66: * application-provided FormJacobian can provide this matrix without interfering with the user's (command-line) choices. */
67: static Mat FieldSplitSchurPre(PC_FieldSplit *jac)
68: {
69: switch (jac->schurpre) {
70: case PC_FIELDSPLIT_SCHUR_PRE_SELF: return jac->schur;
71: case PC_FIELDSPLIT_SCHUR_PRE_SELFP: return jac->schurp;
72: case PC_FIELDSPLIT_SCHUR_PRE_A11: return jac->pmat[1];
73: case PC_FIELDSPLIT_SCHUR_PRE_FULL: /* We calculate this and store it in schur_user */
74: case PC_FIELDSPLIT_SCHUR_PRE_USER: /* Use a user-provided matrix if it is given, otherwise diagonal block */
75: default:
76: return jac->schur_user ? jac->schur_user : jac->pmat[1];
77: }
78: }
81: #include <petscdraw.h>
82: static PetscErrorCode PCView_FieldSplit(PC pc,PetscViewer viewer)
83: {
84: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
85: PetscErrorCode ierr;
86: PetscBool iascii,isdraw;
87: PetscInt i,j;
88: PC_FieldSplitLink ilink = jac->head;
91: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
92: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
93: if (iascii) {
94: if (jac->bs > 0) {
95: PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs);
96: } else {
97: PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D\n",PCCompositeTypes[jac->type],jac->nsplits);
98: }
99: if (pc->useAmat) {
100: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");
101: }
102: if (jac->diag_use_amat) {
103: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for diagonal blocks\n");
104: }
105: if (jac->offdiag_use_amat) {
106: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for off-diagonal blocks\n");
107: }
108: PetscViewerASCIIPrintf(viewer," Solver info for each split is in the following KSP objects:\n");
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: }
129: if (isdraw) {
130: PetscDraw draw;
131: PetscReal x,y,w,wd;
133: PetscViewerDrawGetDraw(viewer,0,&draw);
134: PetscDrawGetCurrentPoint(draw,&x,&y);
135: w = 2*PetscMin(1.0 - x,x);
136: wd = w/(jac->nsplits + 1);
137: x = x - wd*(jac->nsplits-1)/2.0;
138: for (i=0; i<jac->nsplits; i++) {
139: PetscDrawPushCurrentPoint(draw,x,y);
140: KSPView(ilink->ksp,viewer);
141: PetscDrawPopCurrentPoint(draw);
142: x += wd;
143: ilink = ilink->next;
144: }
145: }
146: return(0);
147: }
149: static PetscErrorCode PCView_FieldSplit_Schur(PC pc,PetscViewer viewer)
150: {
151: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
152: PetscErrorCode ierr;
153: PetscBool iascii,isdraw;
154: PetscInt i,j;
155: PC_FieldSplitLink ilink = jac->head;
156: MatSchurComplementAinvType atype;
159: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
160: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
161: if (iascii) {
162: if (jac->bs > 0) {
163: PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, blocksize = %D, factorization %s\n",jac->bs,PCFieldSplitSchurFactTypes[jac->schurfactorization]);
164: } else {
165: PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, factorization %s\n",PCFieldSplitSchurFactTypes[jac->schurfactorization]);
166: }
167: if (pc->useAmat) {
168: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");
169: }
170: switch (jac->schurpre) {
171: case PC_FIELDSPLIT_SCHUR_PRE_SELF:
172: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from S itself\n");
173: break;
174: case PC_FIELDSPLIT_SCHUR_PRE_SELFP:
175: MatSchurComplementGetAinvType(jac->schur,&atype);
176: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from Sp, an assembled approximation to S, which uses A00's %sdiagonal's inverse\n",atype == MAT_SCHUR_COMPLEMENT_AINV_DIAG ? "" : (atype == MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG ? "block " : "lumped "));break;
177: case PC_FIELDSPLIT_SCHUR_PRE_A11:
178: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n");
179: break;
180: case PC_FIELDSPLIT_SCHUR_PRE_FULL:
181: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from the exact Schur complement\n");
182: break;
183: case PC_FIELDSPLIT_SCHUR_PRE_USER:
184: if (jac->schur_user) {
185: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from user provided matrix\n");
186: } else {
187: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n");
188: }
189: break;
190: default:
191: SETERRQ1(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Invalid Schur preconditioning type: %d", jac->schurpre);
192: }
193: PetscViewerASCIIPrintf(viewer," Split info:\n");
194: PetscViewerASCIIPushTab(viewer);
195: for (i=0; i<jac->nsplits; i++) {
196: if (ilink->fields) {
197: PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);
198: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
199: for (j=0; j<ilink->nfields; j++) {
200: if (j > 0) {
201: PetscViewerASCIIPrintf(viewer,",");
202: }
203: PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);
204: }
205: PetscViewerASCIIPrintf(viewer,"\n");
206: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
207: } else {
208: PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);
209: }
210: ilink = ilink->next;
211: }
212: PetscViewerASCIIPrintf(viewer,"KSP solver for A00 block\n");
213: PetscViewerASCIIPushTab(viewer);
214: if (jac->head) {
215: KSPView(jac->head->ksp,viewer);
216: } else {PetscViewerASCIIPrintf(viewer," not yet available\n");}
217: PetscViewerASCIIPopTab(viewer);
218: if (jac->head && jac->kspupper != jac->head->ksp) {
219: PetscViewerASCIIPrintf(viewer,"KSP solver for upper A00 in upper triangular factor \n");
220: PetscViewerASCIIPushTab(viewer);
221: if (jac->kspupper) {KSPView(jac->kspupper,viewer);}
222: else {PetscViewerASCIIPrintf(viewer," not yet available\n");}
223: PetscViewerASCIIPopTab(viewer);
224: }
225: PetscViewerASCIIPrintf(viewer,"KSP solver for S = A11 - A10 inv(A00) A01 \n");
226: PetscViewerASCIIPushTab(viewer);
227: if (jac->kspschur) {
228: KSPView(jac->kspschur,viewer);
229: } else {
230: PetscViewerASCIIPrintf(viewer," not yet available\n");
231: }
232: PetscViewerASCIIPopTab(viewer);
233: PetscViewerASCIIPopTab(viewer);
234: } else if (isdraw && jac->head) {
235: PetscDraw draw;
236: PetscReal x,y,w,wd,h;
237: PetscInt cnt = 2;
238: char str[32];
240: PetscViewerDrawGetDraw(viewer,0,&draw);
241: PetscDrawGetCurrentPoint(draw,&x,&y);
242: if (jac->kspupper != jac->head->ksp) cnt++;
243: w = 2*PetscMin(1.0 - x,x);
244: wd = w/(cnt + 1);
246: PetscSNPrintf(str,32,"Schur fact. %s",PCFieldSplitSchurFactTypes[jac->schurfactorization]);
247: PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
248: y -= h;
249: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_USER && !jac->schur_user) {
250: PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[PC_FIELDSPLIT_SCHUR_PRE_A11]);
251: } else {
252: PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[jac->schurpre]);
253: }
254: PetscDrawStringBoxed(draw,x+wd*(cnt-1)/2.0,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
255: y -= h;
256: x = x - wd*(cnt-1)/2.0;
258: PetscDrawPushCurrentPoint(draw,x,y);
259: KSPView(jac->head->ksp,viewer);
260: PetscDrawPopCurrentPoint(draw);
261: if (jac->kspupper != jac->head->ksp) {
262: x += wd;
263: PetscDrawPushCurrentPoint(draw,x,y);
264: KSPView(jac->kspupper,viewer);
265: PetscDrawPopCurrentPoint(draw);
266: }
267: x += wd;
268: PetscDrawPushCurrentPoint(draw,x,y);
269: KSPView(jac->kspschur,viewer);
270: PetscDrawPopCurrentPoint(draw);
271: }
272: return(0);
273: }
275: /* Precondition: jac->bs is set to a meaningful value */
276: static PetscErrorCode PCFieldSplitSetRuntimeSplits_Private(PC pc)
277: {
279: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
280: PetscInt i,nfields,*ifields,nfields_col,*ifields_col;
281: PetscBool flg,flg_col;
282: char optionname[128],splitname[8],optionname_col[128];
285: PetscMalloc1(jac->bs,&ifields);
286: PetscMalloc1(jac->bs,&ifields_col);
287: for (i=0,flg=PETSC_TRUE;; i++) {
288: PetscSNPrintf(splitname,sizeof(splitname),"%D",i);
289: PetscSNPrintf(optionname,sizeof(optionname),"-pc_fieldsplit_%D_fields",i);
290: PetscSNPrintf(optionname_col,sizeof(optionname_col),"-pc_fieldsplit_%D_fields_col",i);
291: nfields = jac->bs;
292: nfields_col = jac->bs;
293: PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,optionname,ifields,&nfields,&flg);
294: PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,optionname_col,ifields_col,&nfields_col,&flg_col);
295: if (!flg) break;
296: else if (flg && !flg_col) {
297: if (!nfields) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields");
298: PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields);
299: } else {
300: if (!nfields || !nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields");
301: if (nfields != nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Number of row and column fields must match");
302: PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields_col);
303: }
304: }
305: if (i > 0) {
306: /* Makes command-line setting of splits take precedence over setting them in code.
307: Otherwise subsequent calls to PCFieldSplitSetIS() or PCFieldSplitSetFields() would
308: create new splits, which would probably not be what the user wanted. */
309: jac->splitdefined = PETSC_TRUE;
310: }
311: PetscFree(ifields);
312: PetscFree(ifields_col);
313: return(0);
314: }
316: static PetscErrorCode PCFieldSplitSetDefaults(PC pc)
317: {
318: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
319: PetscErrorCode ierr;
320: PC_FieldSplitLink ilink = jac->head;
321: PetscBool fieldsplit_default = PETSC_FALSE,coupling = PETSC_FALSE;
322: PetscInt i;
325: /*
326: Kinda messy, but at least this now uses DMCreateFieldDecomposition().
327: Should probably be rewritten.
328: */
329: if (!ilink) {
330: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_coupling",&coupling,NULL);
331: if (pc->dm && jac->dm_splits && !jac->detect && !coupling) {
332: PetscInt numFields, f, i, j;
333: char **fieldNames;
334: IS *fields;
335: DM *dms;
336: DM subdm[128];
337: PetscBool flg;
339: DMCreateFieldDecomposition(pc->dm, &numFields, &fieldNames, &fields, &dms);
340: /* Allow the user to prescribe the splits */
341: for (i = 0, flg = PETSC_TRUE;; i++) {
342: PetscInt ifields[128];
343: IS compField;
344: char optionname[128], splitname[8];
345: PetscInt nfields = numFields;
347: PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%D_fields", i);
348: PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix, optionname, ifields, &nfields, &flg);
349: if (!flg) break;
350: if (numFields > 128) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot currently support %d > 128 fields", numFields);
351: DMCreateSubDM(pc->dm, nfields, ifields, &compField, &subdm[i]);
352: if (nfields == 1) {
353: PCFieldSplitSetIS(pc, fieldNames[ifields[0]], compField);
354: } else {
355: PetscSNPrintf(splitname, sizeof(splitname), "%D", i);
356: PCFieldSplitSetIS(pc, splitname, compField);
357: }
358: ISDestroy(&compField);
359: for (j = 0; j < nfields; ++j) {
360: f = ifields[j];
361: PetscFree(fieldNames[f]);
362: ISDestroy(&fields[f]);
363: }
364: }
365: if (i == 0) {
366: for (f = 0; f < numFields; ++f) {
367: PCFieldSplitSetIS(pc, fieldNames[f], fields[f]);
368: PetscFree(fieldNames[f]);
369: ISDestroy(&fields[f]);
370: }
371: } else {
372: for (j=0; j<numFields; j++) {
373: DMDestroy(dms+j);
374: }
375: PetscFree(dms);
376: PetscMalloc1(i, &dms);
377: for (j = 0; j < i; ++j) dms[j] = subdm[j];
378: }
379: PetscFree(fieldNames);
380: PetscFree(fields);
381: if (dms) {
382: PetscInfo(pc, "Setting up physics based fieldsplit preconditioner using the embedded DM\n");
383: for (ilink = jac->head, i = 0; ilink; ilink = ilink->next, ++i) {
384: const char *prefix;
385: PetscObjectGetOptionsPrefix((PetscObject)(ilink->ksp),&prefix);
386: PetscObjectSetOptionsPrefix((PetscObject)(dms[i]), prefix);
387: KSPSetDM(ilink->ksp, dms[i]);
388: KSPSetDMActive(ilink->ksp, PETSC_FALSE);
389: PetscObjectIncrementTabLevel((PetscObject)dms[i],(PetscObject)ilink->ksp,0);
390: DMDestroy(&dms[i]);
391: }
392: PetscFree(dms);
393: }
394: } else {
395: if (jac->bs <= 0) {
396: if (pc->pmat) {
397: MatGetBlockSize(pc->pmat,&jac->bs);
398: } else jac->bs = 1;
399: }
401: if (jac->detect) {
402: IS zerodiags,rest;
403: PetscInt nmin,nmax;
405: MatGetOwnershipRange(pc->mat,&nmin,&nmax);
406: MatFindZeroDiagonals(pc->mat,&zerodiags);
407: ISComplement(zerodiags,nmin,nmax,&rest);
408: PCFieldSplitSetIS(pc,"0",rest);
409: PCFieldSplitSetIS(pc,"1",zerodiags);
410: ISDestroy(&zerodiags);
411: ISDestroy(&rest);
412: } else if (coupling) {
413: IS coupling,rest;
414: PetscInt nmin,nmax;
416: MatGetOwnershipRange(pc->mat,&nmin,&nmax);
417: MatFindOffBlockDiagonalEntries(pc->mat,&coupling);
418: ISCreateStride(PetscObjectComm((PetscObject)pc->mat),nmax-nmin,nmin,1,&rest);
419: ISSetIdentity(rest);
420: PCFieldSplitSetIS(pc,"0",rest);
421: PCFieldSplitSetIS(pc,"1",coupling);
422: ISDestroy(&coupling);
423: ISDestroy(&rest);
424: } else {
425: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&fieldsplit_default,NULL);
426: if (!fieldsplit_default) {
427: /* Allow user to set fields from command line, if bs was known at the time of PCSetFromOptions_FieldSplit()
428: then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */
429: PCFieldSplitSetRuntimeSplits_Private(pc);
430: if (jac->splitdefined) {PetscInfo(pc,"Splits defined using the options database\n");}
431: }
432: if ((fieldsplit_default || !jac->splitdefined) && !jac->isrestrict) {
433: PetscInfo(pc,"Using default splitting of fields\n");
434: for (i=0; i<jac->bs; i++) {
435: char splitname[8];
436: PetscSNPrintf(splitname,sizeof(splitname),"%D",i);
437: PCFieldSplitSetFields(pc,splitname,1,&i,&i);
438: }
439: jac->defaultsplit = PETSC_TRUE;
440: }
441: }
442: }
443: } else if (jac->nsplits == 1) {
444: if (ilink->is) {
445: IS is2;
446: PetscInt nmin,nmax;
448: MatGetOwnershipRange(pc->mat,&nmin,&nmax);
449: ISComplement(ilink->is,nmin,nmax,&is2);
450: PCFieldSplitSetIS(pc,"1",is2);
451: ISDestroy(&is2);
452: } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Must provide at least two sets of fields to PCFieldSplit()");
453: }
455: if (jac->nsplits < 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Unhandled case, must have at least two fields, not %d", jac->nsplits);
456: return(0);
457: }
459: PETSC_EXTERN PetscErrorCode PetscOptionsFindPairPrefix_Private(PetscOptions,const char pre[], const char name[],const char *value[],PetscBool *flg);
461: static PetscErrorCode PCSetUp_FieldSplit(PC pc)
462: {
463: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
464: PetscErrorCode ierr;
465: PC_FieldSplitLink ilink;
466: PetscInt i,nsplit;
467: PetscBool sorted, sorted_col;
470: pc->failedreason = PC_NOERROR;
471: PCFieldSplitSetDefaults(pc);
472: nsplit = jac->nsplits;
473: ilink = jac->head;
475: /* get the matrices for each split */
476: if (!jac->issetup) {
477: PetscInt rstart,rend,nslots,bs;
479: jac->issetup = PETSC_TRUE;
481: /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */
482: if (jac->defaultsplit || !ilink->is) {
483: if (jac->bs <= 0) jac->bs = nsplit;
484: }
485: bs = jac->bs;
486: MatGetOwnershipRange(pc->pmat,&rstart,&rend);
487: nslots = (rend - rstart)/bs;
488: for (i=0; i<nsplit; i++) {
489: if (jac->defaultsplit) {
490: ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+i,nsplit,&ilink->is);
491: ISDuplicate(ilink->is,&ilink->is_col);
492: } else if (!ilink->is) {
493: if (ilink->nfields > 1) {
494: PetscInt *ii,*jj,j,k,nfields = ilink->nfields,*fields = ilink->fields,*fields_col = ilink->fields_col;
495: PetscMalloc1(ilink->nfields*nslots,&ii);
496: PetscMalloc1(ilink->nfields*nslots,&jj);
497: for (j=0; j<nslots; j++) {
498: for (k=0; k<nfields; k++) {
499: ii[nfields*j + k] = rstart + bs*j + fields[k];
500: jj[nfields*j + k] = rstart + bs*j + fields_col[k];
501: }
502: }
503: ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,ii,PETSC_OWN_POINTER,&ilink->is);
504: ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,jj,PETSC_OWN_POINTER,&ilink->is_col);
505: ISSetBlockSize(ilink->is, nfields);
506: ISSetBlockSize(ilink->is_col, nfields);
507: } else {
508: ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields[0],bs,&ilink->is);
509: ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields_col[0],bs,&ilink->is_col);
510: }
511: }
512: ISSorted(ilink->is,&sorted);
513: if (ilink->is_col) { ISSorted(ilink->is_col,&sorted_col); }
514: if (!sorted || !sorted_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Fields must be sorted when creating split");
515: ilink = ilink->next;
516: }
517: }
519: ilink = jac->head;
520: if (!jac->pmat) {
521: Vec xtmp;
523: MatCreateVecs(pc->pmat,&xtmp,NULL);
524: PetscMalloc1(nsplit,&jac->pmat);
525: PetscMalloc2(nsplit,&jac->x,nsplit,&jac->y);
526: for (i=0; i<nsplit; i++) {
527: MatNullSpace sp;
529: /* Check for preconditioning matrix attached to IS */
530: PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &jac->pmat[i]);
531: if (jac->pmat[i]) {
532: PetscObjectReference((PetscObject) jac->pmat[i]);
533: if (jac->type == PC_COMPOSITE_SCHUR) {
534: jac->schur_user = jac->pmat[i];
536: PetscObjectReference((PetscObject) jac->schur_user);
537: }
538: } else {
539: const char *prefix;
540: MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->pmat[i]);
541: KSPGetOptionsPrefix(ilink->ksp,&prefix);
542: MatSetOptionsPrefix(jac->pmat[i],prefix);
543: MatViewFromOptions(jac->pmat[i],NULL,"-mat_view");
544: }
545: /* create work vectors for each split */
546: MatCreateVecs(jac->pmat[i],&jac->x[i],&jac->y[i]);
547: ilink->x = jac->x[i]; ilink->y = jac->y[i]; ilink->z = NULL;
548: /* compute scatter contexts needed by multiplicative versions and non-default splits */
549: VecScatterCreate(xtmp,ilink->is,jac->x[i],NULL,&ilink->sctx);
550: PetscObjectQuery((PetscObject) ilink->is, "nearnullspace", (PetscObject*) &sp);
551: if (sp) {
552: MatSetNearNullSpace(jac->pmat[i], sp);
553: }
554: ilink = ilink->next;
555: }
556: VecDestroy(&xtmp);
557: } else {
558: MatReuse scall;
559: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
560: for (i=0; i<nsplit; i++) {
561: MatDestroy(&jac->pmat[i]);
562: }
563: scall = MAT_INITIAL_MATRIX;
564: } else scall = MAT_REUSE_MATRIX;
566: for (i=0; i<nsplit; i++) {
567: Mat pmat;
569: /* Check for preconditioning matrix attached to IS */
570: PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &pmat);
571: if (!pmat) {
572: MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is_col,scall,&jac->pmat[i]);
573: }
574: ilink = ilink->next;
575: }
576: }
577: if (jac->diag_use_amat) {
578: ilink = jac->head;
579: if (!jac->mat) {
580: PetscMalloc1(nsplit,&jac->mat);
581: for (i=0; i<nsplit; i++) {
582: MatCreateSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->mat[i]);
583: ilink = ilink->next;
584: }
585: } else {
586: MatReuse scall;
587: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
588: for (i=0; i<nsplit; i++) {
589: MatDestroy(&jac->mat[i]);
590: }
591: scall = MAT_INITIAL_MATRIX;
592: } else scall = MAT_REUSE_MATRIX;
594: for (i=0; i<nsplit; i++) {
595: if (jac->mat[i]) {MatCreateSubMatrix(pc->mat,ilink->is,ilink->is_col,scall,&jac->mat[i]);}
596: ilink = ilink->next;
597: }
598: }
599: } else {
600: jac->mat = jac->pmat;
601: }
603: /* Check for null space attached to IS */
604: ilink = jac->head;
605: for (i=0; i<nsplit; i++) {
606: MatNullSpace sp;
608: PetscObjectQuery((PetscObject) ilink->is, "nullspace", (PetscObject*) &sp);
609: if (sp) {
610: MatSetNullSpace(jac->mat[i], sp);
611: }
612: ilink = ilink->next;
613: }
615: if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR) {
616: /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */
617: /* FIXME: Can/should we reuse jac->mat whenever (jac->diag_use_amat) is true? */
618: ilink = jac->head;
619: if (nsplit == 2 && jac->type == PC_COMPOSITE_MULTIPLICATIVE) {
620: /* 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 */
621: if (!jac->Afield) {
622: PetscCalloc1(nsplit,&jac->Afield);
623: if (jac->offdiag_use_amat) {
624: MatCreateSubMatrix(pc->mat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);
625: } else {
626: MatCreateSubMatrix(pc->pmat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);
627: }
628: } else {
629: MatReuse scall;
630: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
631: for (i=0; i<nsplit; i++) {
632: MatDestroy(&jac->Afield[1]);
633: }
634: scall = MAT_INITIAL_MATRIX;
635: } else scall = MAT_REUSE_MATRIX;
637: if (jac->offdiag_use_amat) {
638: MatCreateSubMatrix(pc->mat,ilink->next->is,ilink->is,scall,&jac->Afield[1]);
639: } else {
640: MatCreateSubMatrix(pc->pmat,ilink->next->is,ilink->is,scall,&jac->Afield[1]);
641: }
642: }
643: } else {
644: if (!jac->Afield) {
645: PetscMalloc1(nsplit,&jac->Afield);
646: for (i=0; i<nsplit; i++) {
647: if (jac->offdiag_use_amat) {
648: MatCreateSubMatrix(pc->mat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);
649: } else {
650: MatCreateSubMatrix(pc->pmat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);
651: }
652: ilink = ilink->next;
653: }
654: } else {
655: MatReuse scall;
656: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
657: for (i=0; i<nsplit; i++) {
658: MatDestroy(&jac->Afield[i]);
659: }
660: scall = MAT_INITIAL_MATRIX;
661: } else scall = MAT_REUSE_MATRIX;
663: for (i=0; i<nsplit; i++) {
664: if (jac->offdiag_use_amat) {
665: MatCreateSubMatrix(pc->mat,ilink->is,NULL,scall,&jac->Afield[i]);
666: } else {
667: MatCreateSubMatrix(pc->pmat,ilink->is,NULL,scall,&jac->Afield[i]);
668: }
669: ilink = ilink->next;
670: }
671: }
672: }
673: }
675: if (jac->type == PC_COMPOSITE_SCHUR) {
676: IS ccis;
677: PetscBool isspd;
678: PetscInt rstart,rend;
679: char lscname[256];
680: PetscObject LSC_L;
682: if (nsplit != 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"To use Schur complement preconditioner you must have exactly 2 fields");
684: /* If pc->mat is SPD, don't scale by -1 the Schur complement */
685: if (jac->schurscale == (PetscScalar)-1.0) {
686: MatGetOption(pc->pmat,MAT_SPD,&isspd);
687: jac->schurscale = (isspd == PETSC_TRUE) ? 1.0 : -1.0;
688: }
690: /* When extracting off-diagonal submatrices, we take complements from this range */
691: MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);
693: /* need to handle case when one is resetting up the preconditioner */
694: if (jac->schur) {
695: KSP kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper;
697: MatSchurComplementGetKSP(jac->schur, &kspInner);
698: ilink = jac->head;
699: ISComplement(ilink->is_col,rstart,rend,&ccis);
700: if (jac->offdiag_use_amat) {
701: MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);
702: } else {
703: MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);
704: }
705: ISDestroy(&ccis);
706: ilink = ilink->next;
707: ISComplement(ilink->is_col,rstart,rend,&ccis);
708: if (jac->offdiag_use_amat) {
709: MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);
710: } else {
711: MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);
712: }
713: ISDestroy(&ccis);
714: MatSchurComplementUpdateSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);
715: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
716: MatDestroy(&jac->schurp);
717: MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);
718: }
719: if (kspA != kspInner) {
720: KSPSetOperators(kspA,jac->mat[0],jac->pmat[0]);
721: }
722: if (kspUpper != kspA) {
723: KSPSetOperators(kspUpper,jac->mat[0],jac->pmat[0]);
724: }
725: KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));
726: } else {
727: const char *Dprefix;
728: char schurprefix[256], schurmatprefix[256];
729: char schurtestoption[256];
730: MatNullSpace sp;
731: PetscBool flg;
732: KSP kspt;
734: /* extract the A01 and A10 matrices */
735: ilink = jac->head;
736: ISComplement(ilink->is_col,rstart,rend,&ccis);
737: if (jac->offdiag_use_amat) {
738: MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
739: } else {
740: MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
741: }
742: ISDestroy(&ccis);
743: ilink = ilink->next;
744: ISComplement(ilink->is_col,rstart,rend,&ccis);
745: if (jac->offdiag_use_amat) {
746: MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
747: } else {
748: MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
749: }
750: ISDestroy(&ccis);
752: /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */
753: MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur);
754: MatSetType(jac->schur,MATSCHURCOMPLEMENT);
755: MatSchurComplementSetSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);
756: PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
757: MatSetOptionsPrefix(jac->schur,schurmatprefix);
758: MatSchurComplementGetKSP(jac->schur,&kspt);
759: KSPSetOptionsPrefix(kspt,schurmatprefix);
761: /* Note: this is not true in general */
762: MatGetNullSpace(jac->mat[1], &sp);
763: if (sp) {
764: MatSetNullSpace(jac->schur, sp);
765: }
767: PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);
768: PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
769: if (flg) {
770: DM dmInner;
771: KSP kspInner;
772: PC pcInner;
774: MatSchurComplementGetKSP(jac->schur, &kspInner);
775: KSPReset(kspInner);
776: KSPSetOperators(kspInner,jac->mat[0],jac->pmat[0]);
777: PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
778: /* Indent this deeper to emphasize the "inner" nature of this solver. */
779: PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject) pc, 2);
780: PetscObjectIncrementTabLevel((PetscObject)kspInner->pc, (PetscObject) pc, 2);
781: KSPSetOptionsPrefix(kspInner, schurprefix);
783: /* Set DM for new solver */
784: KSPGetDM(jac->head->ksp, &dmInner);
785: KSPSetDM(kspInner, dmInner);
786: KSPSetDMActive(kspInner, PETSC_FALSE);
788: /* Defaults to PCKSP as preconditioner */
789: KSPGetPC(kspInner, &pcInner);
790: PCSetType(pcInner, PCKSP);
791: PCKSPSetKSP(pcInner, jac->head->ksp);
792: } else {
793: /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or
794: * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact,
795: * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for
796: * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make
797: * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used
798: * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */
799: KSPSetType(jac->head->ksp,KSPGMRES);
800: MatSchurComplementSetKSP(jac->schur,jac->head->ksp);
801: }
802: KSPSetOperators(jac->head->ksp,jac->mat[0],jac->pmat[0]);
803: KSPSetFromOptions(jac->head->ksp);
804: MatSetFromOptions(jac->schur);
806: PetscObjectTypeCompare((PetscObject)jac->schur, MATSCHURCOMPLEMENT, &flg);
807: if (flg) { /* Need to do this otherwise PCSetUp_KSP will overwrite the amat of jac->head->ksp */
808: KSP kspInner;
809: PC pcInner;
811: MatSchurComplementGetKSP(jac->schur, &kspInner);
812: KSPGetPC(kspInner, &pcInner);
813: PetscObjectTypeCompare((PetscObject)pcInner, PCKSP, &flg);
814: if (flg) {
815: KSP ksp;
817: PCKSPGetKSP(pcInner, &ksp);
818: if (ksp == jac->head->ksp) {
819: PCSetUseAmat(pcInner, PETSC_TRUE);
820: }
821: }
822: }
823: PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);
824: PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
825: if (flg) {
826: DM dmInner;
828: PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
829: KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper);
830: KSPSetErrorIfNotConverged(jac->kspupper,pc->erroriffailure);
831: KSPSetOptionsPrefix(jac->kspupper, schurprefix);
832: PetscObjectIncrementTabLevel((PetscObject)jac->kspupper, (PetscObject) pc, 1);
833: PetscObjectIncrementTabLevel((PetscObject)jac->kspupper->pc, (PetscObject) pc, 1);
834: KSPGetDM(jac->head->ksp, &dmInner);
835: KSPSetDM(jac->kspupper, dmInner);
836: KSPSetDMActive(jac->kspupper, PETSC_FALSE);
837: KSPSetFromOptions(jac->kspupper);
838: KSPSetOperators(jac->kspupper,jac->mat[0],jac->pmat[0]);
839: VecDuplicate(jac->head->x, &jac->head->z);
840: } else {
841: jac->kspupper = jac->head->ksp;
842: PetscObjectReference((PetscObject) jac->head->ksp);
843: }
845: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
846: MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);
847: }
848: KSPCreate(PetscObjectComm((PetscObject)pc),&jac->kspschur);
849: KSPSetErrorIfNotConverged(jac->kspschur,pc->erroriffailure);
850: PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);
851: PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);
852: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) {
853: PC pcschur;
854: KSPGetPC(jac->kspschur,&pcschur);
855: PCSetType(pcschur,PCNONE);
856: /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */
857: } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) {
858: MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user);
859: }
860: KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));
861: KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix);
862: KSPSetOptionsPrefix(jac->kspschur, Dprefix);
863: /* propagate DM */
864: {
865: DM sdm;
866: KSPGetDM(jac->head->next->ksp, &sdm);
867: if (sdm) {
868: KSPSetDM(jac->kspschur, sdm);
869: KSPSetDMActive(jac->kspschur, PETSC_FALSE);
870: }
871: }
872: /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
873: /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */
874: KSPSetFromOptions(jac->kspschur);
875: }
876: MatAssemblyBegin(jac->schur,MAT_FINAL_ASSEMBLY);
877: MatAssemblyEnd(jac->schur,MAT_FINAL_ASSEMBLY);
879: /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */
880: PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_L",ilink->splitname);
881: PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);
882: if (!LSC_L) {PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);}
883: if (LSC_L) {PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L);}
884: PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_Lp",ilink->splitname);
885: PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);
886: if (!LSC_L) {PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);}
887: if (LSC_L) {PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L);}
888: } else {
889: /* set up the individual splits' PCs */
890: i = 0;
891: ilink = jac->head;
892: while (ilink) {
893: KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i]);
894: /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
895: if (!jac->suboptionsset) {KSPSetFromOptions(ilink->ksp);}
896: i++;
897: ilink = ilink->next;
898: }
899: }
901: jac->suboptionsset = PETSC_TRUE;
902: return(0);
903: }
905: #define FieldSplitSplitSolveAdd(ilink,xx,yy) \
906: (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
907: VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
908: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\
909: KSPSolve(ilink->ksp,ilink->x,ilink->y) || \
910: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\
911: VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \
912: VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE))
914: static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y)
915: {
916: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
917: PetscErrorCode ierr;
918: PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
919: KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;
920: KSPConvergedReason reason;
923: switch (jac->schurfactorization) {
924: case PC_FIELDSPLIT_SCHUR_FACT_DIAG:
925: /* [A00 0; 0 -S], positive definite, suitable for MINRES */
926: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
927: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
928: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
929: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
930: KSPSolve(kspA,ilinkA->x,ilinkA->y);
931: KSPGetConvergedReason(kspA,&reason);
932: if(reason < 0) pc->failedreason = PC_SUBPC_ERROR;
933: PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
934: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
935: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
936: PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
937: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
938: KSPGetConvergedReason(jac->kspschur,&reason);
939: if(reason < 0) pc->failedreason = PC_SUBPC_ERROR;
940: PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
941: VecScale(ilinkD->y,jac->schurscale);
942: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
943: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
944: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
945: break;
946: case PC_FIELDSPLIT_SCHUR_FACT_LOWER:
947: /* [A00 0; A10 S], suitable for left preconditioning */
948: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
949: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
950: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
951: KSPSolve(kspA,ilinkA->x,ilinkA->y);
952: KSPGetConvergedReason(kspA,&reason);
953: if(reason < 0) pc->failedreason = PC_SUBPC_ERROR;
954: PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
955: MatMult(jac->C,ilinkA->y,ilinkD->x);
956: VecScale(ilinkD->x,-1.);
957: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
958: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
959: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
960: PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
961: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
962: KSPGetConvergedReason(jac->kspschur,&reason);
963: if(reason < 0) pc->failedreason = PC_SUBPC_ERROR;
964: PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
965: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
966: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
967: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
968: break;
969: case PC_FIELDSPLIT_SCHUR_FACT_UPPER:
970: /* [A00 A01; 0 S], suitable for right preconditioning */
971: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
972: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
973: PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
974: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
975: KSPGetConvergedReason(jac->kspschur,&reason);
976: if(reason < 0) pc->failedreason = PC_SUBPC_ERROR;
977: PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL); MatMult(jac->B,ilinkD->y,ilinkA->x);
978: VecScale(ilinkA->x,-1.);
979: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);
980: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
981: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);
982: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
983: KSPSolve(kspA,ilinkA->x,ilinkA->y);
984: KSPGetConvergedReason(kspA,&reason);
985: if(reason < 0) pc->failedreason = PC_SUBPC_ERROR;
986: PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
987: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
988: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
989: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
990: break;
991: case PC_FIELDSPLIT_SCHUR_FACT_FULL:
992: /* [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 */
993: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
994: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
995: PetscLogEventBegin(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);
996: KSPSolve(kspLower,ilinkA->x,ilinkA->y);
997: KSPGetConvergedReason(kspLower,&reason);
998: if(reason < 0) pc->failedreason = PC_SUBPC_ERROR;
999: PetscLogEventEnd(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);
1000: MatMult(jac->C,ilinkA->y,ilinkD->x);
1001: VecScale(ilinkD->x,-1.0);
1002: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
1003: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
1005: PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1006: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
1007: KSPGetConvergedReason(jac->kspschur,&reason);
1008: if(reason < 0) pc->failedreason = PC_SUBPC_ERROR;
1009: PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1010: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1011: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1013: if (kspUpper == kspA) {
1014: MatMult(jac->B,ilinkD->y,ilinkA->y);
1015: VecAXPY(ilinkA->x,-1.0,ilinkA->y);
1016: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1017: KSPSolve(kspA,ilinkA->x,ilinkA->y);
1018: KSPGetConvergedReason(kspA,&reason);
1019: if(reason < 0) pc->failedreason = PC_SUBPC_ERROR;
1020: PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1021: } else {
1022: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1023: KSPSolve(kspA,ilinkA->x,ilinkA->y);
1024: PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1025: MatMult(jac->B,ilinkD->y,ilinkA->x);
1026: PetscLogEventBegin(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);
1027: KSPSolve(kspUpper,ilinkA->x,ilinkA->z);
1028: KSPGetConvergedReason(kspUpper,&reason);
1029: if(reason < 0) pc->failedreason = PC_SUBPC_ERROR;
1030: PetscLogEventEnd(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);
1031: VecAXPY(ilinkA->y,-1.0,ilinkA->z);
1032: }
1033: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1034: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1035: }
1036: return(0);
1037: }
1039: static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y)
1040: {
1041: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1042: PetscErrorCode ierr;
1043: PC_FieldSplitLink ilink = jac->head;
1044: PetscInt cnt,bs;
1045: KSPConvergedReason reason;
1048: if (jac->type == PC_COMPOSITE_ADDITIVE) {
1049: if (jac->defaultsplit) {
1050: VecGetBlockSize(x,&bs);
1051: 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);
1052: VecGetBlockSize(y,&bs);
1053: 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);
1054: VecStrideGatherAll(x,jac->x,INSERT_VALUES);
1055: while (ilink) {
1056: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1057: KSPSolve(ilink->ksp,ilink->x,ilink->y);
1058: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1059: KSPGetConvergedReason(ilink->ksp,&reason);
1060: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1061: pc->failedreason = PC_SUBPC_ERROR;
1062: }
1063: ilink = ilink->next;
1064: }
1065: VecStrideScatterAll(jac->y,y,INSERT_VALUES);
1066: } else {
1067: VecSet(y,0.0);
1068: while (ilink) {
1069: FieldSplitSplitSolveAdd(ilink,x,y);
1070: KSPGetConvergedReason(ilink->ksp,&reason);
1071: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1072: pc->failedreason = PC_SUBPC_ERROR;
1073: }
1074: ilink = ilink->next;
1075: }
1076: }
1077: } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) {
1078: VecSet(y,0.0);
1079: /* solve on first block for first block variables */
1080: VecScatterBegin(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);
1081: VecScatterEnd(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);
1082: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1083: KSPSolve(ilink->ksp,ilink->x,ilink->y);
1084: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1085: KSPGetConvergedReason(ilink->ksp,&reason);
1086: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1087: pc->failedreason = PC_SUBPC_ERROR;
1088: }
1089: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1090: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1092: /* compute the residual only onto second block variables using first block variables */
1093: MatMult(jac->Afield[1],ilink->y,ilink->next->x);
1094: ilink = ilink->next;
1095: VecScale(ilink->x,-1.0);
1096: VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1097: VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1099: /* solve on second block variables */
1100: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1101: KSPSolve(ilink->ksp,ilink->x,ilink->y);
1102: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1103: KSPGetConvergedReason(ilink->ksp,&reason);
1104: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1105: pc->failedreason = PC_SUBPC_ERROR;
1106: }
1107: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1108: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1109: } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1110: if (!jac->w1) {
1111: VecDuplicate(x,&jac->w1);
1112: VecDuplicate(x,&jac->w2);
1113: }
1114: VecSet(y,0.0);
1115: FieldSplitSplitSolveAdd(ilink,x,y);
1116: KSPGetConvergedReason(ilink->ksp,&reason);
1117: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1118: pc->failedreason = PC_SUBPC_ERROR;
1119: }
1120: cnt = 1;
1121: while (ilink->next) {
1122: ilink = ilink->next;
1123: /* compute the residual only over the part of the vector needed */
1124: MatMult(jac->Afield[cnt++],y,ilink->x);
1125: VecScale(ilink->x,-1.0);
1126: VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1127: VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1128: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1129: KSPSolve(ilink->ksp,ilink->x,ilink->y);
1130: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1131: KSPGetConvergedReason(ilink->ksp,&reason);
1132: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1133: pc->failedreason = PC_SUBPC_ERROR;
1134: }
1135: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1136: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1137: }
1138: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1139: cnt -= 2;
1140: while (ilink->previous) {
1141: ilink = ilink->previous;
1142: /* compute the residual only over the part of the vector needed */
1143: MatMult(jac->Afield[cnt--],y,ilink->x);
1144: VecScale(ilink->x,-1.0);
1145: VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1146: VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1147: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1148: KSPSolve(ilink->ksp,ilink->x,ilink->y);
1149: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1150: KSPGetConvergedReason(ilink->ksp,&reason);
1151: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1152: pc->failedreason = PC_SUBPC_ERROR;
1153: }
1154: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1155: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1156: }
1157: }
1158: } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type);
1159: return(0);
1160: }
1162: #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \
1163: (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
1164: VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
1165: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \
1166: KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \
1167: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \
1168: VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \
1169: VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE))
1171: static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y)
1172: {
1173: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1174: PetscErrorCode ierr;
1175: PC_FieldSplitLink ilink = jac->head;
1176: PetscInt bs;
1177: KSPConvergedReason reason;
1180: if (jac->type == PC_COMPOSITE_ADDITIVE) {
1181: if (jac->defaultsplit) {
1182: VecGetBlockSize(x,&bs);
1183: 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);
1184: VecGetBlockSize(y,&bs);
1185: 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);
1186: VecStrideGatherAll(x,jac->x,INSERT_VALUES);
1187: while (ilink) {
1188: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1189: KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);
1190: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1191: KSPGetConvergedReason(ilink->ksp,&reason);
1192: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1193: pc->failedreason = PC_SUBPC_ERROR;
1194: }
1195: ilink = ilink->next;
1196: }
1197: VecStrideScatterAll(jac->y,y,INSERT_VALUES);
1198: } else {
1199: VecSet(y,0.0);
1200: while (ilink) {
1201: FieldSplitSplitSolveAddTranspose(ilink,x,y);
1202: KSPGetConvergedReason(ilink->ksp,&reason);
1203: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1204: pc->failedreason = PC_SUBPC_ERROR;
1205: }
1206: ilink = ilink->next;
1207: }
1208: }
1209: } else {
1210: if (!jac->w1) {
1211: VecDuplicate(x,&jac->w1);
1212: VecDuplicate(x,&jac->w2);
1213: }
1214: VecSet(y,0.0);
1215: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1216: FieldSplitSplitSolveAddTranspose(ilink,x,y);
1217: KSPGetConvergedReason(ilink->ksp,&reason);
1218: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1219: pc->failedreason = PC_SUBPC_ERROR;
1220: }
1221: while (ilink->next) {
1222: ilink = ilink->next;
1223: MatMultTranspose(pc->mat,y,jac->w1);
1224: VecWAXPY(jac->w2,-1.0,jac->w1,x);
1225: FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1226: }
1227: while (ilink->previous) {
1228: ilink = ilink->previous;
1229: MatMultTranspose(pc->mat,y,jac->w1);
1230: VecWAXPY(jac->w2,-1.0,jac->w1,x);
1231: FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1232: }
1233: } else {
1234: while (ilink->next) { /* get to last entry in linked list */
1235: ilink = ilink->next;
1236: }
1237: FieldSplitSplitSolveAddTranspose(ilink,x,y);
1238: KSPGetConvergedReason(ilink->ksp,&reason);
1239: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1240: pc->failedreason = PC_SUBPC_ERROR;
1241: }
1242: while (ilink->previous) {
1243: ilink = ilink->previous;
1244: MatMultTranspose(pc->mat,y,jac->w1);
1245: VecWAXPY(jac->w2,-1.0,jac->w1,x);
1246: FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1247: }
1248: }
1249: }
1250: return(0);
1251: }
1253: static PetscErrorCode PCReset_FieldSplit(PC pc)
1254: {
1255: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1256: PetscErrorCode ierr;
1257: PC_FieldSplitLink ilink = jac->head,next;
1260: while (ilink) {
1261: KSPDestroy(&ilink->ksp);
1262: VecDestroy(&ilink->x);
1263: VecDestroy(&ilink->y);
1264: VecDestroy(&ilink->z);
1265: VecScatterDestroy(&ilink->sctx);
1266: ISDestroy(&ilink->is);
1267: ISDestroy(&ilink->is_col);
1268: PetscFree(ilink->splitname);
1269: PetscFree(ilink->fields);
1270: PetscFree(ilink->fields_col);
1271: next = ilink->next;
1272: PetscFree(ilink);
1273: ilink = next;
1274: }
1275: jac->head = NULL;
1276: PetscFree2(jac->x,jac->y);
1277: if (jac->mat && jac->mat != jac->pmat) {
1278: MatDestroyMatrices(jac->nsplits,&jac->mat);
1279: } else if (jac->mat) {
1280: jac->mat = NULL;
1281: }
1282: if (jac->pmat) {MatDestroyMatrices(jac->nsplits,&jac->pmat);}
1283: if (jac->Afield) {MatDestroyMatrices(jac->nsplits,&jac->Afield);}
1284: jac->nsplits = 0;
1285: VecDestroy(&jac->w1);
1286: VecDestroy(&jac->w2);
1287: MatDestroy(&jac->schur);
1288: MatDestroy(&jac->schurp);
1289: MatDestroy(&jac->schur_user);
1290: KSPDestroy(&jac->kspschur);
1291: KSPDestroy(&jac->kspupper);
1292: MatDestroy(&jac->B);
1293: MatDestroy(&jac->C);
1294: jac->isrestrict = PETSC_FALSE;
1295: return(0);
1296: }
1298: static PetscErrorCode PCDestroy_FieldSplit(PC pc)
1299: {
1300: PetscErrorCode ierr;
1303: PCReset_FieldSplit(pc);
1304: PetscFree(pc->data);
1305: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurGetSubKSP_C",NULL);
1306: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL);
1307: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL);
1308: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL);
1309: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL);
1310: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL);
1311: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL);
1312: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL);
1313: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL);
1314: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",NULL);
1315: return(0);
1316: }
1318: static PetscErrorCode PCSetFromOptions_FieldSplit(PetscOptionItems *PetscOptionsObject,PC pc)
1319: {
1320: PetscErrorCode ierr;
1321: PetscInt bs;
1322: PetscBool flg;
1323: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1324: PCCompositeType ctype;
1327: PetscOptionsHead(PetscOptionsObject,"FieldSplit options");
1328: PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL);
1329: PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);
1330: if (flg) {
1331: PCFieldSplitSetBlockSize(pc,bs);
1332: }
1333: jac->diag_use_amat = pc->useAmat;
1334: 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);
1335: jac->offdiag_use_amat = pc->useAmat;
1336: 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);
1337: PetscOptionsBool("-pc_fieldsplit_detect_saddle_point","Form 2-way split by detecting zero diagonal entries", "PCFieldSplitSetDetectSaddlePoint",jac->detect,&jac->detect,NULL);
1338: PCFieldSplitSetDetectSaddlePoint(pc,jac->detect); /* Sets split type and Schur PC type */
1339: PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);
1340: if (flg) {
1341: PCFieldSplitSetType(pc,ctype);
1342: }
1343: /* Only setup fields once */
1344: if ((jac->bs > 0) && (jac->nsplits == 0)) {
1345: /* only allow user to set fields from command line if bs is already known.
1346: otherwise user can set them in PCFieldSplitSetDefaults() */
1347: PCFieldSplitSetRuntimeSplits_Private(pc);
1348: if (jac->splitdefined) {PetscInfo(pc,"Splits defined using the options database\n");}
1349: }
1350: if (jac->type == PC_COMPOSITE_SCHUR) {
1351: PetscOptionsGetEnum(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);
1352: if (flg) {PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");}
1353: 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);
1354: PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSetSchurPre",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL);
1355: PetscOptionsScalar("-pc_fieldsplit_schur_scale","Scale Schur complement","PCFieldSplitSetSchurScale",jac->schurscale,&jac->schurscale,NULL);
1356: }
1357: PetscOptionsTail();
1358: return(0);
1359: }
1361: /*------------------------------------------------------------------------------------*/
1363: static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1364: {
1365: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1366: PetscErrorCode ierr;
1367: PC_FieldSplitLink ilink,next = jac->head;
1368: char prefix[128];
1369: PetscInt i;
1372: if (jac->splitdefined) {
1373: PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1374: return(0);
1375: }
1376: for (i=0; i<n; i++) {
1377: if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs);
1378: if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]);
1379: }
1380: PetscNew(&ilink);
1381: if (splitname) {
1382: PetscStrallocpy(splitname,&ilink->splitname);
1383: } else {
1384: PetscMalloc1(3,&ilink->splitname);
1385: PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);
1386: }
1387: 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 */
1388: PetscMalloc1(n,&ilink->fields);
1389: PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));
1390: PetscMalloc1(n,&ilink->fields_col);
1391: PetscMemcpy(ilink->fields_col,fields_col,n*sizeof(PetscInt));
1393: ilink->nfields = n;
1394: ilink->next = NULL;
1395: KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1396: KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);
1397: PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1398: KSPSetType(ilink->ksp,KSPPREONLY);
1399: PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);
1401: PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);
1402: KSPSetOptionsPrefix(ilink->ksp,prefix);
1404: if (!next) {
1405: jac->head = ilink;
1406: ilink->previous = NULL;
1407: } else {
1408: while (next->next) {
1409: next = next->next;
1410: }
1411: next->next = ilink;
1412: ilink->previous = next;
1413: }
1414: jac->nsplits++;
1415: return(0);
1416: }
1418: static PetscErrorCode PCFieldSplitSchurGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp)
1419: {
1420: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1424: *subksp = NULL;
1425: if (n) *n = 0;
1426: if (jac->type == PC_COMPOSITE_SCHUR) {
1427: PetscInt nn;
1429: if (!jac->schur) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitSchurGetSubKSP()");
1430: if (jac->nsplits != 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Unexpected number of splits %D != 2",jac->nsplits);
1431: nn = jac->nsplits + (jac->kspupper != jac->head->ksp ? 1 : 0);
1432: PetscMalloc1(nn,subksp);
1433: (*subksp)[0] = jac->head->ksp;
1434: (*subksp)[1] = jac->kspschur;
1435: if (jac->kspupper != jac->head->ksp) (*subksp)[2] = jac->kspupper;
1436: if (n) *n = nn;
1437: }
1438: return(0);
1439: }
1441: static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp)
1442: {
1443: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1447: if (!jac->schur) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitGetSubKSP()");
1448: PetscMalloc1(jac->nsplits,subksp);
1449: MatSchurComplementGetKSP(jac->schur,*subksp);
1451: (*subksp)[1] = jac->kspschur;
1452: if (n) *n = jac->nsplits;
1453: return(0);
1454: }
1456: static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp)
1457: {
1458: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1459: PetscErrorCode ierr;
1460: PetscInt cnt = 0;
1461: PC_FieldSplitLink ilink = jac->head;
1464: PetscMalloc1(jac->nsplits,subksp);
1465: while (ilink) {
1466: (*subksp)[cnt++] = ilink->ksp;
1467: ilink = ilink->next;
1468: }
1469: 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);
1470: if (n) *n = jac->nsplits;
1471: return(0);
1472: }
1474: /*@C
1475: PCFieldSplitRestrictIS - Restricts the fieldsplit ISs to be within a given IS.
1477: Input Parameters:
1478: + pc - the preconditioner context
1479: + is - the index set that defines the indices to which the fieldsplit is to be restricted
1481: Level: advanced
1483: @*/
1484: PetscErrorCode PCFieldSplitRestrictIS(PC pc,IS isy)
1485: {
1491: PetscTryMethod(pc,"PCFieldSplitRestrictIS_C",(PC,IS),(pc,isy));
1492: return(0);
1493: }
1496: static PetscErrorCode PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy)
1497: {
1498: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1499: PetscErrorCode ierr;
1500: PC_FieldSplitLink ilink = jac->head, next;
1501: PetscInt localsize,size,sizez,i;
1502: const PetscInt *ind, *indz;
1503: PetscInt *indc, *indcz;
1504: PetscBool flg;
1507: ISGetLocalSize(isy,&localsize);
1508: MPI_Scan(&localsize,&size,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)isy));
1509: size -= localsize;
1510: while(ilink) {
1511: IS isrl,isr;
1512: PC subpc;
1513: ISEmbed(ilink->is, isy, PETSC_TRUE, &isrl);
1514: ISGetLocalSize(isrl,&localsize);
1515: PetscMalloc1(localsize,&indc);
1516: ISGetIndices(isrl,&ind);
1517: PetscMemcpy(indc,ind,localsize*sizeof(PetscInt));
1518: ISRestoreIndices(isrl,&ind);
1519: ISDestroy(&isrl);
1520: for (i=0; i<localsize; i++) *(indc+i) += size;
1521: ISCreateGeneral(PetscObjectComm((PetscObject)isy),localsize,indc,PETSC_OWN_POINTER,&isr);
1522: PetscObjectReference((PetscObject)isr);
1523: ISDestroy(&ilink->is);
1524: ilink->is = isr;
1525: PetscObjectReference((PetscObject)isr);
1526: ISDestroy(&ilink->is_col);
1527: ilink->is_col = isr;
1528: ISDestroy(&isr);
1529: KSPGetPC(ilink->ksp, &subpc);
1530: PetscObjectTypeCompare((PetscObject)subpc,PCFIELDSPLIT,&flg);
1531: if(flg) {
1532: IS iszl,isz;
1533: MPI_Comm comm;
1534: ISGetLocalSize(ilink->is,&localsize);
1535: comm = PetscObjectComm((PetscObject)ilink->is);
1536: ISEmbed(isy, ilink->is, PETSC_TRUE, &iszl);
1537: MPI_Scan(&localsize,&sizez,1,MPIU_INT,MPI_SUM,comm);
1538: sizez -= localsize;
1539: ISGetLocalSize(iszl,&localsize);
1540: PetscMalloc1(localsize,&indcz);
1541: ISGetIndices(iszl,&indz);
1542: PetscMemcpy(indcz,indz,localsize*sizeof(PetscInt));
1543: ISRestoreIndices(iszl,&indz);
1544: ISDestroy(&iszl);
1545: for (i=0; i<localsize; i++) *(indcz+i) += sizez;
1546: ISCreateGeneral(comm,localsize,indcz,PETSC_OWN_POINTER,&isz);
1547: PCFieldSplitRestrictIS(subpc,isz);
1548: ISDestroy(&isz);
1549: }
1550: next = ilink->next;
1551: ilink = next;
1552: }
1553: jac->isrestrict = PETSC_TRUE;
1554: return(0);
1555: }
1557: static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is)
1558: {
1559: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1560: PetscErrorCode ierr;
1561: PC_FieldSplitLink ilink, next = jac->head;
1562: char prefix[128];
1565: if (jac->splitdefined) {
1566: PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1567: return(0);
1568: }
1569: PetscNew(&ilink);
1570: if (splitname) {
1571: PetscStrallocpy(splitname,&ilink->splitname);
1572: } else {
1573: PetscMalloc1(8,&ilink->splitname);
1574: PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);
1575: }
1576: 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 */
1577: PetscObjectReference((PetscObject)is);
1578: ISDestroy(&ilink->is);
1579: ilink->is = is;
1580: PetscObjectReference((PetscObject)is);
1581: ISDestroy(&ilink->is_col);
1582: ilink->is_col = is;
1583: ilink->next = NULL;
1584: KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1585: KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);
1586: PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1587: KSPSetType(ilink->ksp,KSPPREONLY);
1588: PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);
1590: PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);
1591: KSPSetOptionsPrefix(ilink->ksp,prefix);
1593: if (!next) {
1594: jac->head = ilink;
1595: ilink->previous = NULL;
1596: } else {
1597: while (next->next) {
1598: next = next->next;
1599: }
1600: next->next = ilink;
1601: ilink->previous = next;
1602: }
1603: jac->nsplits++;
1604: return(0);
1605: }
1607: /*@C
1608: PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner
1610: Logically Collective on PC
1612: Input Parameters:
1613: + pc - the preconditioner context
1614: . splitname - name of this split, if NULL the number of the split is used
1615: . n - the number of fields in this split
1616: - fields - the fields in this split
1618: Level: intermediate
1620: Notes:
1621: Use PCFieldSplitSetIS() to set a completely general set of indices as a field.
1623: The PCFieldSplitSetFields() is for defining fields as strided blocks. For example, if the block
1624: 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
1625: 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x....
1626: where the numbered entries indicate what is in the field.
1628: This function is called once per split (it creates a new split each time). Solve options
1629: for this split will be available under the prefix -fieldsplit_SPLITNAME_.
1631: Developer Note: This routine does not actually create the IS representing the split, that is delayed
1632: until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be
1633: available when this routine is called.
1635: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS()
1637: @*/
1638: PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1639: {
1645: if (n < 1) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname);
1647: PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt*,const PetscInt*),(pc,splitname,n,fields,fields_col));
1648: return(0);
1649: }
1651: /*@
1652: PCFieldSplitSetDiagUseAmat - set flag indicating whether to extract diagonal blocks from Amat (rather than Pmat)
1654: Logically Collective on PC
1656: Input Parameters:
1657: + pc - the preconditioner object
1658: - flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from
1660: Options Database:
1661: . -pc_fieldsplit_diag_use_amat
1663: Level: intermediate
1665: .seealso: PCFieldSplitGetDiagUseAmat(), PCFieldSplitSetOffDiagUseAmat(), PCFIELDSPLIT
1667: @*/
1668: PetscErrorCode PCFieldSplitSetDiagUseAmat(PC pc,PetscBool flg)
1669: {
1670: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1671: PetscBool isfs;
1676: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1677: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1678: jac->diag_use_amat = flg;
1679: return(0);
1680: }
1682: /*@
1683: PCFieldSplitGetDiagUseAmat - get the flag indicating whether to extract diagonal blocks from Amat (rather than Pmat)
1685: Logically Collective on PC
1687: Input Parameters:
1688: . pc - the preconditioner object
1690: Output Parameters:
1691: . flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from
1694: Level: intermediate
1696: .seealso: PCFieldSplitSetDiagUseAmat(), PCFieldSplitGetOffDiagUseAmat(), PCFIELDSPLIT
1698: @*/
1699: PetscErrorCode PCFieldSplitGetDiagUseAmat(PC pc,PetscBool *flg)
1700: {
1701: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1702: PetscBool isfs;
1708: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1709: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1710: *flg = jac->diag_use_amat;
1711: return(0);
1712: }
1714: /*@
1715: PCFieldSplitSetOffDiagUseAmat - set flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat)
1717: Logically Collective on PC
1719: Input Parameters:
1720: + pc - the preconditioner object
1721: - flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from
1723: Options Database:
1724: . -pc_fieldsplit_off_diag_use_amat
1726: Level: intermediate
1728: .seealso: PCFieldSplitGetOffDiagUseAmat(), PCFieldSplitSetDiagUseAmat(), PCFIELDSPLIT
1730: @*/
1731: PetscErrorCode PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg)
1732: {
1733: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1734: PetscBool isfs;
1739: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1740: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1741: jac->offdiag_use_amat = flg;
1742: return(0);
1743: }
1745: /*@
1746: PCFieldSplitGetOffDiagUseAmat - get the flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat)
1748: Logically Collective on PC
1750: Input Parameters:
1751: . pc - the preconditioner object
1753: Output Parameters:
1754: . flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from
1757: Level: intermediate
1759: .seealso: PCFieldSplitSetOffDiagUseAmat(), PCFieldSplitGetDiagUseAmat(), PCFIELDSPLIT
1761: @*/
1762: PetscErrorCode PCFieldSplitGetOffDiagUseAmat(PC pc,PetscBool *flg)
1763: {
1764: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1765: PetscBool isfs;
1771: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1772: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1773: *flg = jac->offdiag_use_amat;
1774: return(0);
1775: }
1779: /*@C
1780: PCFieldSplitSetIS - Sets the exact elements for field
1782: Logically Collective on PC
1784: Input Parameters:
1785: + pc - the preconditioner context
1786: . splitname - name of this split, if NULL the number of the split is used
1787: - is - the index set that defines the vector elements in this field
1790: Notes:
1791: Use PCFieldSplitSetFields(), for fields defined by strided types.
1793: This function is called once per split (it creates a new split each time). Solve options
1794: for this split will be available under the prefix -fieldsplit_SPLITNAME_.
1796: Level: intermediate
1798: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize()
1800: @*/
1801: PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is)
1802: {
1809: PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));
1810: return(0);
1811: }
1813: /*@C
1814: PCFieldSplitGetIS - Retrieves the elements for a field as an IS
1816: Logically Collective on PC
1818: Input Parameters:
1819: + pc - the preconditioner context
1820: - splitname - name of this split
1822: Output Parameter:
1823: - is - the index set that defines the vector elements in this field, or NULL if the field is not found
1825: Level: intermediate
1827: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS()
1829: @*/
1830: PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is)
1831: {
1838: {
1839: PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;
1840: PC_FieldSplitLink ilink = jac->head;
1841: PetscBool found;
1843: *is = NULL;
1844: while (ilink) {
1845: PetscStrcmp(ilink->splitname, splitname, &found);
1846: if (found) {
1847: *is = ilink->is;
1848: break;
1849: }
1850: ilink = ilink->next;
1851: }
1852: }
1853: return(0);
1854: }
1856: /*@
1857: PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the
1858: fieldsplit preconditioner. If not set the matrix block size is used.
1860: Logically Collective on PC
1862: Input Parameters:
1863: + pc - the preconditioner context
1864: - bs - the block size
1866: Level: intermediate
1868: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields()
1870: @*/
1871: PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs)
1872: {
1878: PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));
1879: return(0);
1880: }
1882: /*@C
1883: PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits
1885: Collective on KSP
1887: Input Parameter:
1888: . pc - the preconditioner context
1890: Output Parameters:
1891: + n - the number of splits
1892: - subksp - the array of KSP contexts
1894: Note:
1895: After PCFieldSplitGetSubKSP() the array of KSPs is to be freed by the user with PetscFree()
1896: (not the KSP just the array that contains them).
1898: You must call PCSetUp() before calling PCFieldSplitGetSubKSP().
1900: If the fieldsplit is of type PC_COMPOSITE_SCHUR, it returns the KSP object used inside the
1901: Schur complement and the KSP object used to iterate over the Schur complement.
1902: To access all the KSP objects used in PC_COMPOSITE_SCHUR, use PCFieldSplitSchurGetSubKSP()
1904: Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs.
1905: You can call PCFieldSplitGetSubKSP(pc,n,PETSC_NULL_KSP,ierr) to determine how large the
1906: KSP array must be.
1909: Level: advanced
1911: .seealso: PCFIELDSPLIT
1912: @*/
1913: PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
1914: {
1920: PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));
1921: return(0);
1922: }
1924: /*@C
1925: PCFieldSplitSchurGetSubKSP - Gets the KSP contexts used inside the Schur complement based PCFIELDSPLIT
1927: Collective on KSP
1929: Input Parameter:
1930: . pc - the preconditioner context
1932: Output Parameters:
1933: + n - the number of splits
1934: - subksp - the array of KSP contexts
1936: Note:
1937: After PCFieldSplitSchurGetSubKSP() the array of KSPs is to be freed by the user with PetscFree()
1938: (not the KSP just the array that contains them).
1940: You must call PCSetUp() before calling PCFieldSplitSchurGetSubKSP().
1942: If the fieldsplit type is of type PC_COMPOSITE_SCHUR, it returns (in order)
1943: - the KSP used for the (1,1) block
1944: - the KSP used for the Schur complement (not the one used for the interior Schur solver)
1945: - the KSP used for the (1,1) block in the upper triangular factor (if different from that of the (1,1) block).
1947: It returns a null array if the fieldsplit is not of type PC_COMPOSITE_SCHUR; in this case, you should use PCFieldSplitGetSubKSP().
1949: Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs.
1950: You can call PCFieldSplitSchurGetSubKSP(pc,n,PETSC_NULL_KSP,ierr) to determine how large the
1951: KSP array must be.
1953: Level: advanced
1955: .seealso: PCFIELDSPLIT
1956: @*/
1957: PetscErrorCode PCFieldSplitSchurGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
1958: {
1964: PetscUseMethod(pc,"PCFieldSplitSchurGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));
1965: return(0);
1966: }
1968: /*@
1969: PCFieldSplitSetSchurPre - Indicates what operator is used to construct the preconditioner for the Schur complement.
1970: A11 matrix. Otherwise no preconditioner is used.
1972: Collective on PC
1974: Input Parameters:
1975: + pc - the preconditioner context
1976: . 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
1977: PC_FIELDSPLIT_SCHUR_PRE_SELFP, and PC_FIELDSPLIT_SCHUR_PRE_FULL
1978: - userpre - matrix to use for preconditioning, or NULL
1980: Options Database:
1981: . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11. See notes for meaning of various arguments
1983: Notes:
1984: $ If ptype is
1985: $ a11 then the preconditioner for the Schur complement is generated from the block diagonal part of the preconditioner
1986: $ matrix associated with the Schur complement (i.e. A11), not the Schur complement matrix
1987: $ self the preconditioner for the Schur complement is generated from the symbolic representation of the Schur complement matrix:
1988: $ The only preconditioner that currently works with this symbolic respresentation matrix object is the PCLSC
1989: $ preconditioner
1990: $ user then the preconditioner for the Schur complement is generated from the user provided matrix (pre argument
1991: $ to this function).
1992: $ selfp then the preconditioning for the Schur complement is generated from an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01
1993: $ This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be
1994: $ lumped before extracting the diagonal using the additional option -fieldsplit_1_mat_schur_complement_ainv_type lump
1995: $ full then the preconditioner for the Schur complement is generated from the exact Schur complement matrix representation computed internally by PCFIELDSPLIT (this is expensive)
1996: $ useful mostly as a test that the Schur complement approach can work for your problem
1998: When solving a saddle point problem, where the A11 block is identically zero, using a11 as the ptype only makes sense
1999: with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and
2000: -fieldsplit_1_pc_type lsc which uses the least squares commutator to compute a preconditioner for the Schur complement.
2002: Level: intermediate
2004: .seealso: PCFieldSplitGetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType,
2005: MatSchurComplementSetAinvType(), PCLSC
2007: @*/
2008: PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
2009: {
2014: PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));
2015: return(0);
2016: }
2018: PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) {return PCFieldSplitSetSchurPre(pc,ptype,pre);} /* Deprecated name */
2020: /*@
2021: PCFieldSplitGetSchurPre - For Schur complement fieldsplit, determine how the Schur complement will be
2022: preconditioned. See PCFieldSplitSetSchurPre() for details.
2024: Logically Collective on PC
2026: Input Parameters:
2027: . pc - the preconditioner context
2029: Output Parameters:
2030: + ptype - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11, PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_PRE_USER
2031: - userpre - matrix to use for preconditioning (with PC_FIELDSPLIT_PRE_USER), or NULL
2033: Level: intermediate
2035: .seealso: PCFieldSplitSetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC
2037: @*/
2038: PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
2039: {
2044: PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre));
2045: return(0);
2046: }
2048: /*@
2049: PCFieldSplitSchurGetS - extract the MatSchurComplement object used by this PC in case it needs to be configured separately
2051: Not collective
2053: Input Parameter:
2054: . pc - the preconditioner context
2056: Output Parameter:
2057: . S - the Schur complement matrix
2059: Notes:
2060: This matrix should not be destroyed using MatDestroy(); rather, use PCFieldSplitSchurRestoreS().
2062: Level: advanced
2064: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurRestoreS()
2066: @*/
2067: PetscErrorCode PCFieldSplitSchurGetS(PC pc,Mat *S)
2068: {
2070: const char* t;
2071: PetscBool isfs;
2072: PC_FieldSplit *jac;
2076: PetscObjectGetType((PetscObject)pc,&t);
2077: PetscStrcmp(t,PCFIELDSPLIT,&isfs);
2078: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
2079: jac = (PC_FieldSplit*)pc->data;
2080: if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
2081: if (S) *S = jac->schur;
2082: return(0);
2083: }
2085: /*@
2086: PCFieldSplitSchurRestoreS - restores the MatSchurComplement object used by this PC
2088: Not collective
2090: Input Parameters:
2091: + pc - the preconditioner context
2092: . S - the Schur complement matrix
2094: Level: advanced
2096: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurGetS()
2098: @*/
2099: PetscErrorCode PCFieldSplitSchurRestoreS(PC pc,Mat *S)
2100: {
2102: const char* t;
2103: PetscBool isfs;
2104: PC_FieldSplit *jac;
2108: PetscObjectGetType((PetscObject)pc,&t);
2109: PetscStrcmp(t,PCFIELDSPLIT,&isfs);
2110: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
2111: jac = (PC_FieldSplit*)pc->data;
2112: if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
2113: if (!S || *S != jac->schur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten");
2114: return(0);
2115: }
2118: static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
2119: {
2120: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2124: jac->schurpre = ptype;
2125: if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) {
2126: MatDestroy(&jac->schur_user);
2127: jac->schur_user = pre;
2128: PetscObjectReference((PetscObject)jac->schur_user);
2129: }
2130: return(0);
2131: }
2133: static PetscErrorCode PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
2134: {
2135: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2138: *ptype = jac->schurpre;
2139: *pre = jac->schur_user;
2140: return(0);
2141: }
2143: /*@
2144: PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain in the preconditioner
2146: Collective on PC
2148: Input Parameters:
2149: + pc - the preconditioner context
2150: - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default
2152: Options Database:
2153: . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full
2156: Level: intermediate
2158: Notes:
2159: The FULL factorization is
2161: $ (A B) = (1 0) (A 0) (1 Ainv*B) = L D U
2162: $ (C E) (C*Ainv 1) (0 S) (0 1 )
2164: where S = E - C*Ainv*B. In practice, the full factorization is applied via block triangular solves with the grouping L*(D*U). UPPER uses D*U, LOWER uses L*D,
2165: and DIAG is the diagonal part with the sign of S flipped (because this makes the preconditioner positive definite for many formulations, thus allowing the use of KSPMINRES). Sign flipping of S can be turned off with PCFieldSplitSetSchurScale().
2167: $ If A and S are solved exactly
2168: $ *) FULL factorization is a direct solver.
2169: $ *) The preconditioned operator with LOWER or UPPER has all eigenvalues equal to 1 and minimal polynomial of degree 2, so KSPGMRES converges in 2 iterations.
2170: $ *) With DIAG, the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so KSPGMRES converges in at most 4 iterations.
2172: If the iteration count is very low, consider using KSPFGMRES or KSPGCR which can use one less preconditioner
2173: application in this case. Note that the preconditioned operator may be highly non-normal, so such fast convergence may not be observed in practice.
2175: For symmetric problems in which A is positive definite and S is negative definite, DIAG can be used with KSPMINRES.
2177: Note that a flexible method like KSPFGMRES or KSPGCR must be used if the fieldsplit preconditioner is nonlinear (e.g. a few iterations of a Krylov method is used to solve with A or S).
2179: References:
2180: + 1. - Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000).
2181: - 2. - Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001).
2183: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCFieldSplitSetSchurScale()
2184: @*/
2185: PetscErrorCode PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype)
2186: {
2191: PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));
2192: return(0);
2193: }
2195: static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype)
2196: {
2197: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2200: jac->schurfactorization = ftype;
2201: return(0);
2202: }
2204: /*@
2205: PCFieldSplitSetSchurScale - Controls the sign flip of S for PC_FIELDSPLIT_SCHUR_FACT_DIAG.
2207: Collective on PC
2209: Input Parameters:
2210: + pc - the preconditioner context
2211: - scale - scaling factor for the Schur complement
2213: Options Database:
2214: . -pc_fieldsplit_schur_scale - default is -1.0
2216: Level: intermediate
2218: .seealso: PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurFactType, PCFieldSplitSetSchurScale()
2219: @*/
2220: PetscErrorCode PCFieldSplitSetSchurScale(PC pc,PetscScalar scale)
2221: {
2227: PetscTryMethod(pc,"PCFieldSplitSetSchurScale_C",(PC,PetscScalar),(pc,scale));
2228: return(0);
2229: }
2231: static PetscErrorCode PCFieldSplitSetSchurScale_FieldSplit(PC pc,PetscScalar scale)
2232: {
2233: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2236: jac->schurscale = scale;
2237: return(0);
2238: }
2240: /*@C
2241: PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement
2243: Collective on KSP
2245: Input Parameter:
2246: . pc - the preconditioner context
2248: Output Parameters:
2249: + A00 - the (0,0) block
2250: . A01 - the (0,1) block
2251: . A10 - the (1,0) block
2252: - A11 - the (1,1) block
2254: Level: advanced
2256: .seealso: PCFIELDSPLIT
2257: @*/
2258: PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11)
2259: {
2260: PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;
2264: if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach.");
2265: if (A00) *A00 = jac->pmat[0];
2266: if (A01) *A01 = jac->B;
2267: if (A10) *A10 = jac->C;
2268: if (A11) *A11 = jac->pmat[1];
2269: return(0);
2270: }
2272: static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type)
2273: {
2274: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2278: jac->type = type;
2279: if (type == PC_COMPOSITE_SCHUR) {
2280: pc->ops->apply = PCApply_FieldSplit_Schur;
2281: pc->ops->view = PCView_FieldSplit_Schur;
2283: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);
2284: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit);
2285: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit);
2286: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);
2287: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",PCFieldSplitSetSchurScale_FieldSplit);
2289: } else {
2290: pc->ops->apply = PCApply_FieldSplit;
2291: pc->ops->view = PCView_FieldSplit;
2293: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2294: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",0);
2295: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",0);
2296: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);
2297: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",0);
2298: }
2299: return(0);
2300: }
2302: static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs)
2303: {
2304: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2307: if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs);
2308: 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);
2309: jac->bs = bs;
2310: return(0);
2311: }
2313: /*@
2314: PCFieldSplitSetType - Sets the type of fieldsplit preconditioner.
2316: Collective on PC
2318: Input Parameter:
2319: . pc - the preconditioner context
2320: . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR
2322: Options Database Key:
2323: . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type
2325: Level: Intermediate
2327: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
2329: .seealso: PCCompositeSetType()
2331: @*/
2332: PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type)
2333: {
2338: PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));
2339: return(0);
2340: }
2342: /*@
2343: PCFieldSplitGetType - Gets the type of fieldsplit preconditioner.
2345: Not collective
2347: Input Parameter:
2348: . pc - the preconditioner context
2350: Output Parameter:
2351: . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR
2353: Level: Intermediate
2355: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
2356: .seealso: PCCompositeSetType()
2357: @*/
2358: PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type)
2359: {
2360: PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;
2365: *type = jac->type;
2366: return(0);
2367: }
2369: /*@
2370: PCFieldSplitSetDMSplits - Flags whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible.
2372: Logically Collective
2374: Input Parameters:
2375: + pc - the preconditioner context
2376: - flg - boolean indicating whether to use field splits defined by the DM
2378: Options Database Key:
2379: . -pc_fieldsplit_dm_splits
2381: Level: Intermediate
2383: .keywords: PC, DM, composite preconditioner, additive, multiplicative
2385: .seealso: PCFieldSplitGetDMSplits()
2387: @*/
2388: PetscErrorCode PCFieldSplitSetDMSplits(PC pc,PetscBool flg)
2389: {
2390: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2391: PetscBool isfs;
2397: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2398: if (isfs) {
2399: jac->dm_splits = flg;
2400: }
2401: return(0);
2402: }
2405: /*@
2406: PCFieldSplitGetDMSplits - Returns flag indicating whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible.
2408: Logically Collective
2410: Input Parameter:
2411: . pc - the preconditioner context
2413: Output Parameter:
2414: . flg - boolean indicating whether to use field splits defined by the DM
2416: Level: Intermediate
2418: .keywords: PC, DM, composite preconditioner, additive, multiplicative
2420: .seealso: PCFieldSplitSetDMSplits()
2422: @*/
2423: PetscErrorCode PCFieldSplitGetDMSplits(PC pc,PetscBool* flg)
2424: {
2425: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2426: PetscBool isfs;
2432: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2433: if (isfs) {
2434: if(flg) *flg = jac->dm_splits;
2435: }
2436: return(0);
2437: }
2439: /*@
2440: PCFieldSplitGetDetectSaddlePoint - Returns flag indicating whether PCFieldSplit will attempt to automatically determine fields based on zero diagonal entries.
2442: Logically Collective
2444: Input Parameter:
2445: . pc - the preconditioner context
2447: Output Parameter:
2448: . flg - boolean indicating whether to detect fields or not
2450: Level: Intermediate
2452: .seealso: PCFIELDSPLIT, PCFieldSplitSetDetectSaddlePoint()
2454: @*/
2455: PetscErrorCode PCFieldSplitGetDetectSaddlePoint(PC pc,PetscBool *flg)
2456: {
2457: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2460: *flg = jac->detect;
2461: return(0);
2462: }
2464: /*@
2465: PCFieldSplitSetDetectSaddlePoint - Sets flag indicating whether PCFieldSplit will attempt to automatically determine fields based on zero diagonal entries.
2467: Logically Collective
2469: Notes:
2470: Also sets the split type to PC_COMPOSITE_SCHUR (see PCFieldSplitSetType()) and the Schur preconditioner type to PC_FIELDSPLIT_SCHUR_PRE_SELF (see PCFieldSplitSetSchurPre()).
2472: Input Parameter:
2473: . pc - the preconditioner context
2475: Output Parameter:
2476: . flg - boolean indicating whether to detect fields or not
2478: Options Database Key:
2479: . -pc_fieldsplit_detect_saddle_point
2481: Level: Intermediate
2483: .seealso: PCFIELDSPLIT, PCFieldSplitSetDetectSaddlePoint(), PCFieldSplitSetType(), PCFieldSplitSetSchurPre()
2485: @*/
2486: PetscErrorCode PCFieldSplitSetDetectSaddlePoint(PC pc,PetscBool flg)
2487: {
2488: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2492: jac->detect = flg;
2493: if (jac->detect) {
2494: PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);
2495: PCFieldSplitSetSchurPre(pc,PC_FIELDSPLIT_SCHUR_PRE_SELF,NULL);
2496: }
2497: return(0);
2498: }
2500: /* -------------------------------------------------------------------------------------*/
2501: /*MC
2502: PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual
2503: fields or groups of fields. See the users manual section "Solving Block Matrices" for more details.
2505: To set options on the solvers for each block append -fieldsplit_ to all the PC
2506: options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1
2508: To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP()
2509: and set the options directly on the resulting KSP object
2511: Level: intermediate
2513: Options Database Keys:
2514: + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split
2515: . -pc_fieldsplit_default - automatically add any fields to additional splits that have not
2516: been supplied explicitly by -pc_fieldsplit_%d_fields
2517: . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields)
2518: . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur> - type of relaxation or factorization splitting
2519: . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11; see PCFieldSplitSetSchurPre()
2520: . -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero diagonal and uses Schur complement with no preconditioner as the solver
2522: - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_
2523: for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields
2525: Notes:
2526: Use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS()
2527: to define a field by an arbitrary collection of entries.
2529: If no fields are set the default is used. The fields are defined by entries strided by bs,
2530: beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(),
2531: if this is not called the block size defaults to the blocksize of the second matrix passed
2532: to KSPSetOperators()/PCSetOperators().
2534: $ For the Schur complement preconditioner if J = ( A00 A01 )
2535: $ ( A10 A11 )
2536: $ the preconditioner using full factorization is
2537: $ ( I -ksp(A00) A01 ) ( inv(A00) 0 ) ( I 0 )
2538: $ ( 0 I ) ( 0 ksp(S) ) ( -A10 ksp(A00) I )
2539: where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. S is the Schur complement
2540: $ S = A11 - A10 ksp(A00) A01
2541: 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
2542: in providing the SECOND split or 1 if not give). For PCFieldSplitGetSubKSP() when field number is 0,
2543: it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default
2544: A11 is used to construct a preconditioner for S, use PCFieldSplitSetSchurPre() for all the possible ways to construct the preconditioner for S.
2546: The factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above,
2547: diag gives
2548: $ ( inv(A00) 0 )
2549: $ ( 0 -ksp(S) )
2550: note that slightly counter intuitively there is a negative in front of the ksp(S) so that the preconditioner is positive definite. For SPD matrices J, the sign flip
2551: can be turned off with PCFieldSplitSetSchurScale() or by command line -pc_fieldsplit_schur_scale 1.0. The lower factorization is the inverse of
2552: $ ( A00 0 )
2553: $ ( A10 S )
2554: where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of
2555: $ ( A00 A01 )
2556: $ ( 0 S )
2557: where again the inverses of A00 and S are applied using KSPs.
2559: If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS
2560: is used automatically for a second block.
2562: The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1.
2563: Generally it should be used with the AIJ format.
2565: The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see,
2566: for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT
2567: inside a smoother resulting in "Distributive Smoothers".
2569: Concepts: physics based preconditioners, block preconditioners
2571: There is a nice discussion of block preconditioners in
2573: [El08] A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier-Stokes equations
2574: Howard Elman, V.E. Howle, John Shadid, Robert Shuttleworth, Ray Tuminaro, Journal of Computational Physics 227 (2008) 1790--1808
2575: http://chess.cs.umd.edu/~elman/papers/tax.pdf
2577: The Constrained Pressure Preconditioner (CPR) can be implemented using PCCOMPOSITE with PCGALERKIN. CPR first solves an R A P subsystem, updates the
2578: residual on all variables (PCCompositeSetType(pc,PC_COMPOSITE_MULTIPLICATIVE)), and then applies a simple ILU like preconditioner on all the variables.
2580: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC,
2581: PCFieldSplitGetSubKSP(), PCFieldSplitSchurGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre(),
2582: MatSchurComplementSetAinvType(), PCFieldSplitSetSchurScale(),
2583: PCFieldSplitSetDetectSaddlePoint()
2584: M*/
2586: PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc)
2587: {
2589: PC_FieldSplit *jac;
2592: PetscNewLog(pc,&jac);
2594: jac->bs = -1;
2595: jac->nsplits = 0;
2596: jac->type = PC_COMPOSITE_MULTIPLICATIVE;
2597: jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */
2598: jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL;
2599: jac->schurscale = -1.0;
2600: jac->dm_splits = PETSC_TRUE;
2601: jac->detect = PETSC_FALSE;
2603: pc->data = (void*)jac;
2605: pc->ops->apply = PCApply_FieldSplit;
2606: pc->ops->applytranspose = PCApplyTranspose_FieldSplit;
2607: pc->ops->setup = PCSetUp_FieldSplit;
2608: pc->ops->reset = PCReset_FieldSplit;
2609: pc->ops->destroy = PCDestroy_FieldSplit;
2610: pc->ops->setfromoptions = PCSetFromOptions_FieldSplit;
2611: pc->ops->view = PCView_FieldSplit;
2612: pc->ops->applyrichardson = 0;
2614: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurGetSubKSP_C",PCFieldSplitSchurGetSubKSP_FieldSplit);
2615: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2616: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);
2617: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);
2618: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);
2619: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);
2620: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",PCFieldSplitRestrictIS_FieldSplit);
2621: return(0);
2622: }