Actual source code: fieldsplit.c
petsc-3.6.1 2015-08-06
2: #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/
3: #include <petscdm.h>
6: const char *const PCFieldSplitSchurPreTypes[] = {"SELF","SELFP","A11","USER","FULL","PCFieldSplitSchurPreType","PC_FIELDSPLIT_SCHUR_PRE_",0};
7: const char *const PCFieldSplitSchurFactTypes[] = {"DIAG","LOWER","UPPER","FULL","PCFieldSplitSchurFactType","PC_FIELDSPLIT_SCHUR_FACT_",0};
9: typedef struct _PC_FieldSplitLink *PC_FieldSplitLink;
10: struct _PC_FieldSplitLink {
11: KSP ksp;
12: Vec x,y,z;
13: char *splitname;
14: PetscInt nfields;
15: PetscInt *fields,*fields_col;
16: VecScatter sctx;
17: IS is,is_col;
18: PC_FieldSplitLink next,previous;
19: };
21: typedef struct {
22: PCCompositeType type;
23: PetscBool defaultsplit; /* Flag for a system with a set of 'k' scalar fields with the same layout (and bs = k) */
24: PetscBool splitdefined; /* Flag is set after the splits have been defined, to prevent more splits from being added */
25: PetscInt bs; /* Block size for IS and Mat structures */
26: PetscInt nsplits; /* Number of field divisions defined */
27: Vec *x,*y,w1,w2;
28: Mat *mat; /* The diagonal block for each split */
29: Mat *pmat; /* The preconditioning diagonal block for each split */
30: Mat *Afield; /* The rows of the matrix associated with each split */
31: PetscBool issetup;
33: /* Only used when Schur complement preconditioning is used */
34: Mat B; /* The (0,1) block */
35: Mat C; /* The (1,0) block */
36: Mat schur; /* The Schur complement S = A11 - A10 A00^{-1} A01, the KSP here, kspinner, is H_1 in [El08] */
37: Mat schurp; /* Assembled approximation to S built by MatSchurComplement to be used as a preconditioning matrix when solving with S */
38: Mat schur_user; /* User-provided preconditioning matrix for the Schur complement */
39: PCFieldSplitSchurPreType schurpre; /* Determines which preconditioning matrix is used for the Schur complement */
40: PCFieldSplitSchurFactType schurfactorization;
41: KSP kspschur; /* The solver for S */
42: KSP kspupper; /* The solver for A in the upper diagonal part of the factorization (H_2 in [El08]) */
43: PC_FieldSplitLink head;
44: PetscBool reset; /* indicates PCReset() has been last called on this object, hack */
45: PetscBool suboptionsset; /* Indicates that the KSPSetFromOptions() has been called on the sub-KSPs */
46: PetscBool dm_splits; /* Whether to use DMCreateFieldDecomposition() whenever possible */
47: PetscBool diag_use_amat; /* Whether to extract diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */
48: PetscBool offdiag_use_amat; /* Whether to extract off-diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */
49: } PC_FieldSplit;
51: /*
52: Notes: there is no particular reason that pmat, x, and y are stored as arrays in PC_FieldSplit instead of
53: inside PC_FieldSplitLink, just historical. If you want to be able to add new fields after already using the
54: PC you could change this.
55: */
57: /* This helper is so that setting a user-provided preconditioning matrix is orthogonal to choosing to use it. This way the
58: * application-provided FormJacobian can provide this matrix without interfering with the user's (command-line) choices. */
59: static Mat FieldSplitSchurPre(PC_FieldSplit *jac)
60: {
61: switch (jac->schurpre) {
62: case PC_FIELDSPLIT_SCHUR_PRE_SELF: return jac->schur;
63: case PC_FIELDSPLIT_SCHUR_PRE_SELFP: return jac->schurp;
64: case PC_FIELDSPLIT_SCHUR_PRE_A11: return jac->pmat[1];
65: case PC_FIELDSPLIT_SCHUR_PRE_FULL: /* We calculate this and store it in schur_user */
66: case PC_FIELDSPLIT_SCHUR_PRE_USER: /* Use a user-provided matrix if it is given, otherwise diagonal block */
67: default:
68: return jac->schur_user ? jac->schur_user : jac->pmat[1];
69: }
70: }
73: #include <petscdraw.h>
76: static PetscErrorCode PCView_FieldSplit(PC pc,PetscViewer viewer)
77: {
78: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
79: PetscErrorCode ierr;
80: PetscBool iascii,isdraw;
81: PetscInt i,j;
82: PC_FieldSplitLink ilink = jac->head;
85: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
86: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
87: if (iascii) {
88: if (jac->bs > 0) {
89: PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs);
90: } else {
91: PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D\n",PCCompositeTypes[jac->type],jac->nsplits);
92: }
93: if (pc->useAmat) {
94: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");
95: }
96: if (jac->diag_use_amat) {
97: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for diagonal blocks\n");
98: }
99: if (jac->offdiag_use_amat) {
100: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for off-diagonal blocks\n");
101: }
102: PetscViewerASCIIPrintf(viewer," Solver info for each split is in the following KSP objects:\n");
103: PetscViewerASCIIPushTab(viewer);
104: for (i=0; i<jac->nsplits; i++) {
105: if (ilink->fields) {
106: PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);
107: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
108: for (j=0; j<ilink->nfields; j++) {
109: if (j > 0) {
110: PetscViewerASCIIPrintf(viewer,",");
111: }
112: PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);
113: }
114: PetscViewerASCIIPrintf(viewer,"\n");
115: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
116: } else {
117: PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);
118: }
119: KSPView(ilink->ksp,viewer);
120: ilink = ilink->next;
121: }
122: PetscViewerASCIIPopTab(viewer);
123: }
125: if (isdraw) {
126: PetscDraw draw;
127: PetscReal x,y,w,wd;
129: PetscViewerDrawGetDraw(viewer,0,&draw);
130: PetscDrawGetCurrentPoint(draw,&x,&y);
131: w = 2*PetscMin(1.0 - x,x);
132: wd = w/(jac->nsplits + 1);
133: x = x - wd*(jac->nsplits-1)/2.0;
134: for (i=0; i<jac->nsplits; i++) {
135: PetscDrawPushCurrentPoint(draw,x,y);
136: KSPView(ilink->ksp,viewer);
137: PetscDrawPopCurrentPoint(draw);
138: x += wd;
139: ilink = ilink->next;
140: }
141: }
142: return(0);
143: }
147: static PetscErrorCode PCView_FieldSplit_Schur(PC pc,PetscViewer viewer)
148: {
149: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
150: PetscErrorCode ierr;
151: PetscBool iascii,isdraw;
152: PetscInt i,j;
153: PC_FieldSplitLink ilink = jac->head;
156: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
157: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
158: if (iascii) {
159: if (jac->bs > 0) {
160: PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, blocksize = %D, factorization %s\n",jac->bs,PCFieldSplitSchurFactTypes[jac->schurfactorization]);
161: } else {
162: PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, factorization %s\n",PCFieldSplitSchurFactTypes[jac->schurfactorization]);
163: }
164: if (pc->useAmat) {
165: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");
166: }
167: switch (jac->schurpre) {
168: case PC_FIELDSPLIT_SCHUR_PRE_SELF:
169: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from S itself\n");break;
170: case PC_FIELDSPLIT_SCHUR_PRE_SELFP:
171: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from Sp, an assembled approximation to S, which uses (lumped, if requested) A00's diagonal's inverse\n");break;
172: case PC_FIELDSPLIT_SCHUR_PRE_A11:
173: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n");break;
174: case PC_FIELDSPLIT_SCHUR_PRE_FULL:
175: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from the exact Schur complement\n");break;
176: case PC_FIELDSPLIT_SCHUR_PRE_USER:
177: if (jac->schur_user) {
178: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from user provided matrix\n");
179: } else {
180: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n");
181: }
182: break;
183: default:
184: SETERRQ1(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Invalid Schur preconditioning type: %d", jac->schurpre);
185: }
186: PetscViewerASCIIPrintf(viewer," Split info:\n");
187: PetscViewerASCIIPushTab(viewer);
188: for (i=0; i<jac->nsplits; i++) {
189: if (ilink->fields) {
190: PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);
191: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
192: for (j=0; j<ilink->nfields; j++) {
193: if (j > 0) {
194: PetscViewerASCIIPrintf(viewer,",");
195: }
196: PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);
197: }
198: PetscViewerASCIIPrintf(viewer,"\n");
199: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
200: } else {
201: PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);
202: }
203: ilink = ilink->next;
204: }
205: PetscViewerASCIIPrintf(viewer,"KSP solver for A00 block\n");
206: PetscViewerASCIIPushTab(viewer);
207: if (jac->head) {
208: KSPView(jac->head->ksp,viewer);
209: } else {PetscViewerASCIIPrintf(viewer," not yet available\n");}
210: PetscViewerASCIIPopTab(viewer);
211: if (jac->head && jac->kspupper != jac->head->ksp) {
212: PetscViewerASCIIPrintf(viewer,"KSP solver for upper A00 in upper triangular factor \n");
213: PetscViewerASCIIPushTab(viewer);
214: if (jac->kspupper) {KSPView(jac->kspupper,viewer);}
215: else {PetscViewerASCIIPrintf(viewer," not yet available\n");}
216: PetscViewerASCIIPopTab(viewer);
217: }
218: PetscViewerASCIIPrintf(viewer,"KSP solver for S = A11 - A10 inv(A00) A01 \n");
219: PetscViewerASCIIPushTab(viewer);
220: if (jac->kspschur) {
221: KSPView(jac->kspschur,viewer);
222: } else {
223: PetscViewerASCIIPrintf(viewer," not yet available\n");
224: }
225: PetscViewerASCIIPopTab(viewer);
226: PetscViewerASCIIPopTab(viewer);
227: } else if (isdraw && jac->head) {
228: PetscDraw draw;
229: PetscReal x,y,w,wd,h;
230: PetscInt cnt = 2;
231: char str[32];
233: PetscViewerDrawGetDraw(viewer,0,&draw);
234: PetscDrawGetCurrentPoint(draw,&x,&y);
235: if (jac->kspupper != jac->head->ksp) cnt++;
236: w = 2*PetscMin(1.0 - x,x);
237: wd = w/(cnt + 1);
239: PetscSNPrintf(str,32,"Schur fact. %s",PCFieldSplitSchurFactTypes[jac->schurfactorization]);
240: PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
241: y -= h;
242: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_USER && !jac->schur_user) {
243: PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[PC_FIELDSPLIT_SCHUR_PRE_A11]);
244: } else {
245: PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[jac->schurpre]);
246: }
247: PetscDrawStringBoxed(draw,x+wd*(cnt-1)/2.0,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
248: y -= h;
249: x = x - wd*(cnt-1)/2.0;
251: PetscDrawPushCurrentPoint(draw,x,y);
252: KSPView(jac->head->ksp,viewer);
253: PetscDrawPopCurrentPoint(draw);
254: if (jac->kspupper != jac->head->ksp) {
255: x += wd;
256: PetscDrawPushCurrentPoint(draw,x,y);
257: KSPView(jac->kspupper,viewer);
258: PetscDrawPopCurrentPoint(draw);
259: }
260: x += wd;
261: PetscDrawPushCurrentPoint(draw,x,y);
262: KSPView(jac->kspschur,viewer);
263: PetscDrawPopCurrentPoint(draw);
264: }
265: return(0);
266: }
270: /* Precondition: jac->bs is set to a meaningful value */
271: static PetscErrorCode PCFieldSplitSetRuntimeSplits_Private(PC pc)
272: {
274: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
275: PetscInt i,nfields,*ifields,nfields_col,*ifields_col;
276: PetscBool flg,flg_col;
277: char optionname[128],splitname[8],optionname_col[128];
280: PetscMalloc1(jac->bs,&ifields);
281: PetscMalloc1(jac->bs,&ifields_col);
282: for (i=0,flg=PETSC_TRUE;; i++) {
283: PetscSNPrintf(splitname,sizeof(splitname),"%D",i);
284: PetscSNPrintf(optionname,sizeof(optionname),"-pc_fieldsplit_%D_fields",i);
285: PetscSNPrintf(optionname_col,sizeof(optionname_col),"-pc_fieldsplit_%D_fields_col",i);
286: nfields = jac->bs;
287: nfields_col = jac->bs;
288: PetscOptionsGetIntArray(((PetscObject)pc)->prefix,optionname,ifields,&nfields,&flg);
289: PetscOptionsGetIntArray(((PetscObject)pc)->prefix,optionname_col,ifields_col,&nfields_col,&flg_col);
290: if (!flg) break;
291: else if (flg && !flg_col) {
292: if (!nfields) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields");
293: PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields);
294: } else {
295: if (!nfields || !nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields");
296: if (nfields != nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Number of row and column fields must match");
297: PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields_col);
298: }
299: }
300: if (i > 0) {
301: /* Makes command-line setting of splits take precedence over setting them in code.
302: Otherwise subsequent calls to PCFieldSplitSetIS() or PCFieldSplitSetFields() would
303: create new splits, which would probably not be what the user wanted. */
304: jac->splitdefined = PETSC_TRUE;
305: }
306: PetscFree(ifields);
307: PetscFree(ifields_col);
308: return(0);
309: }
313: static PetscErrorCode PCFieldSplitSetDefaults(PC pc)
314: {
315: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
316: PetscErrorCode ierr;
317: PC_FieldSplitLink ilink = jac->head;
318: PetscBool fieldsplit_default = PETSC_FALSE,stokes = PETSC_FALSE,coupling = PETSC_FALSE;
319: PetscInt i;
322: /*
323: Kinda messy, but at least this now uses DMCreateFieldDecomposition() even with jac->reset.
324: Should probably be rewritten.
325: */
326: if (!ilink || jac->reset) {
327: PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);
328: PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_coupling",&coupling,NULL);
329: if (pc->dm && jac->dm_splits && !stokes && !coupling) {
330: PetscInt numFields, f, i, j;
331: char **fieldNames;
332: IS *fields;
333: DM *dms;
334: DM subdm[128];
335: PetscBool flg;
337: DMCreateFieldDecomposition(pc->dm, &numFields, &fieldNames, &fields, &dms);
338: /* Allow the user to prescribe the splits */
339: for (i = 0, flg = PETSC_TRUE;; i++) {
340: PetscInt ifields[128];
341: IS compField;
342: char optionname[128], splitname[8];
343: PetscInt nfields = numFields;
345: PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%D_fields", i);
346: PetscOptionsGetIntArray(((PetscObject) pc)->prefix, optionname, ifields, &nfields, &flg);
347: if (!flg) break;
348: if (numFields > 128) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot currently support %d > 128 fields", numFields);
349: DMCreateSubDM(pc->dm, nfields, ifields, &compField, &subdm[i]);
350: if (nfields == 1) {
351: PCFieldSplitSetIS(pc, fieldNames[ifields[0]], compField);
352: /* PetscPrintf(PetscObjectComm((PetscObject)pc), "%s Field Indices:", fieldNames[ifields[0]]);
353: ISView(compField, NULL); */
354: } else {
355: PetscSNPrintf(splitname, sizeof(splitname), "%D", i);
356: PCFieldSplitSetIS(pc, splitname, compField);
357: /* PetscPrintf(PetscObjectComm((PetscObject)pc), "%s Field Indices:", splitname);
358: ISView(compField, NULL); */
359: }
360: ISDestroy(&compField);
361: for (j = 0; j < nfields; ++j) {
362: f = ifields[j];
363: PetscFree(fieldNames[f]);
364: ISDestroy(&fields[f]);
365: }
366: }
367: if (i == 0) {
368: for (f = 0; f < numFields; ++f) {
369: PCFieldSplitSetIS(pc, fieldNames[f], fields[f]);
370: PetscFree(fieldNames[f]);
371: ISDestroy(&fields[f]);
372: }
373: } else {
374: for (j=0; j<numFields; j++) {
375: DMDestroy(dms+j);
376: }
377: PetscFree(dms);
378: PetscMalloc1(i, &dms);
379: for (j = 0; j < i; ++j) dms[j] = subdm[j];
380: }
381: PetscFree(fieldNames);
382: PetscFree(fields);
383: if (dms) {
384: PetscInfo(pc, "Setting up physics based fieldsplit preconditioner using the embedded DM\n");
385: for (ilink = jac->head, i = 0; ilink; ilink = ilink->next, ++i) {
386: const char *prefix;
387: PetscObjectGetOptionsPrefix((PetscObject)(ilink->ksp),&prefix);
388: PetscObjectSetOptionsPrefix((PetscObject)(dms[i]), prefix);
389: KSPSetDM(ilink->ksp, dms[i]);
390: KSPSetDMActive(ilink->ksp, PETSC_FALSE);
391: PetscObjectIncrementTabLevel((PetscObject)dms[i],(PetscObject)ilink->ksp,0);
392: DMDestroy(&dms[i]);
393: }
394: PetscFree(dms);
395: }
396: } else {
397: if (jac->bs <= 0) {
398: if (pc->pmat) {
399: MatGetBlockSize(pc->pmat,&jac->bs);
400: } else jac->bs = 1;
401: }
403: if (stokes) {
404: IS zerodiags,rest;
405: PetscInt nmin,nmax;
407: MatGetOwnershipRange(pc->mat,&nmin,&nmax);
408: MatFindZeroDiagonals(pc->mat,&zerodiags);
409: ISComplement(zerodiags,nmin,nmax,&rest);
410: if (jac->reset) {
411: jac->head->is = rest;
412: jac->head->next->is = zerodiags;
413: } else {
414: PCFieldSplitSetIS(pc,"0",rest);
415: PCFieldSplitSetIS(pc,"1",zerodiags);
416: }
417: ISDestroy(&zerodiags);
418: ISDestroy(&rest);
419: } else if (coupling) {
420: IS coupling,rest;
421: PetscInt nmin,nmax;
423: MatGetOwnershipRange(pc->mat,&nmin,&nmax);
424: MatFindOffBlockDiagonalEntries(pc->mat,&coupling);
425: ISCreateStride(PetscObjectComm((PetscObject)pc->mat),nmax-nmin,nmin,1,&rest);
426: ISSetIdentity(rest);
427: if (jac->reset) {
428: jac->head->is = coupling;
429: jac->head->next->is = rest;
430: } else {
431: PCFieldSplitSetIS(pc,"0",coupling);
432: PCFieldSplitSetIS(pc,"1",rest);
433: }
434: ISDestroy(&coupling);
435: ISDestroy(&rest);
436: } else {
437: if (jac->reset) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cases not yet handled when PCReset() was used");
438: PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&fieldsplit_default,NULL);
439: if (!fieldsplit_default) {
440: /* Allow user to set fields from command line, if bs was known at the time of PCSetFromOptions_FieldSplit()
441: then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */
442: PCFieldSplitSetRuntimeSplits_Private(pc);
443: if (jac->splitdefined) {PetscInfo(pc,"Splits defined using the options database\n");}
444: }
445: if (fieldsplit_default || !jac->splitdefined) {
446: PetscInfo(pc,"Using default splitting of fields\n");
447: for (i=0; i<jac->bs; i++) {
448: char splitname[8];
449: PetscSNPrintf(splitname,sizeof(splitname),"%D",i);
450: PCFieldSplitSetFields(pc,splitname,1,&i,&i);
451: }
452: jac->defaultsplit = PETSC_TRUE;
453: }
454: }
455: }
456: } else if (jac->nsplits == 1) {
457: if (ilink->is) {
458: IS is2;
459: PetscInt nmin,nmax;
461: MatGetOwnershipRange(pc->mat,&nmin,&nmax);
462: ISComplement(ilink->is,nmin,nmax,&is2);
463: PCFieldSplitSetIS(pc,"1",is2);
464: ISDestroy(&is2);
465: } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Must provide at least two sets of fields to PCFieldSplit()");
466: }
469: if (jac->nsplits < 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Unhandled case, must have at least two fields, not %d", jac->nsplits);
470: return(0);
471: }
473: PETSC_EXTERN PetscErrorCode PetscOptionsFindPairPrefix_Private(const char pre[], const char name[], char *value[], PetscBool *flg);
477: static PetscErrorCode PCSetUp_FieldSplit(PC pc)
478: {
479: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
480: PetscErrorCode ierr;
481: PC_FieldSplitLink ilink;
482: PetscInt i,nsplit;
483: PetscBool sorted, sorted_col;
486: PCFieldSplitSetDefaults(pc);
487: nsplit = jac->nsplits;
488: ilink = jac->head;
490: /* get the matrices for each split */
491: if (!jac->issetup) {
492: PetscInt rstart,rend,nslots,bs;
494: jac->issetup = PETSC_TRUE;
496: /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */
497: if (jac->defaultsplit || !ilink->is) {
498: if (jac->bs <= 0) jac->bs = nsplit;
499: }
500: bs = jac->bs;
501: MatGetOwnershipRange(pc->pmat,&rstart,&rend);
502: nslots = (rend - rstart)/bs;
503: for (i=0; i<nsplit; i++) {
504: if (jac->defaultsplit) {
505: ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+i,nsplit,&ilink->is);
506: ISDuplicate(ilink->is,&ilink->is_col);
507: } else if (!ilink->is) {
508: if (ilink->nfields > 1) {
509: PetscInt *ii,*jj,j,k,nfields = ilink->nfields,*fields = ilink->fields,*fields_col = ilink->fields_col;
510: PetscMalloc1(ilink->nfields*nslots,&ii);
511: PetscMalloc1(ilink->nfields*nslots,&jj);
512: for (j=0; j<nslots; j++) {
513: for (k=0; k<nfields; k++) {
514: ii[nfields*j + k] = rstart + bs*j + fields[k];
515: jj[nfields*j + k] = rstart + bs*j + fields_col[k];
516: }
517: }
518: ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,ii,PETSC_OWN_POINTER,&ilink->is);
519: ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,jj,PETSC_OWN_POINTER,&ilink->is_col);
520: ISSetBlockSize(ilink->is, nfields);
521: ISSetBlockSize(ilink->is_col, nfields);
522: } else {
523: ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields[0],bs,&ilink->is);
524: ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields_col[0],bs,&ilink->is_col);
525: }
526: }
527: ISSorted(ilink->is,&sorted);
528: if (ilink->is_col) { ISSorted(ilink->is_col,&sorted_col); }
529: if (!sorted || !sorted_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Fields must be sorted when creating split");
530: ilink = ilink->next;
531: }
532: }
534: ilink = jac->head;
535: if (!jac->pmat) {
536: Vec xtmp;
538: MatCreateVecs(pc->pmat,&xtmp,NULL);
539: PetscMalloc1(nsplit,&jac->pmat);
540: PetscMalloc2(nsplit,&jac->x,nsplit,&jac->y);
541: for (i=0; i<nsplit; i++) {
542: MatNullSpace sp;
544: /* Check for preconditioning matrix attached to IS */
545: PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &jac->pmat[i]);
546: if (jac->pmat[i]) {
547: PetscObjectReference((PetscObject) jac->pmat[i]);
548: if (jac->type == PC_COMPOSITE_SCHUR) {
549: jac->schur_user = jac->pmat[i];
551: PetscObjectReference((PetscObject) jac->schur_user);
552: }
553: } else {
554: const char *prefix;
555: MatGetSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->pmat[i]);
556: KSPGetOptionsPrefix(ilink->ksp,&prefix);
557: MatSetOptionsPrefix(jac->pmat[i],prefix);
558: MatViewFromOptions(jac->pmat[i],NULL,"-mat_view");
559: }
560: /* create work vectors for each split */
561: MatCreateVecs(jac->pmat[i],&jac->x[i],&jac->y[i]);
562: ilink->x = jac->x[i]; ilink->y = jac->y[i]; ilink->z = NULL;
563: /* compute scatter contexts needed by multiplicative versions and non-default splits */
564: VecScatterCreate(xtmp,ilink->is,jac->x[i],NULL,&ilink->sctx);
565: PetscObjectQuery((PetscObject) ilink->is, "nearnullspace", (PetscObject*) &sp);
566: if (sp) {
567: MatSetNearNullSpace(jac->pmat[i], sp);
568: }
569: ilink = ilink->next;
570: }
571: VecDestroy(&xtmp);
572: } else {
573: for (i=0; i<nsplit; i++) {
574: Mat pmat;
576: /* Check for preconditioning matrix attached to IS */
577: PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &pmat);
578: if (!pmat) {
579: MatGetSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_REUSE_MATRIX,&jac->pmat[i]);
580: }
581: ilink = ilink->next;
582: }
583: }
584: if (jac->diag_use_amat) {
585: ilink = jac->head;
586: if (!jac->mat) {
587: PetscMalloc1(nsplit,&jac->mat);
588: for (i=0; i<nsplit; i++) {
589: MatGetSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->mat[i]);
590: ilink = ilink->next;
591: }
592: } else {
593: for (i=0; i<nsplit; i++) {
594: if (jac->mat[i]) {MatGetSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_REUSE_MATRIX,&jac->mat[i]);}
595: ilink = ilink->next;
596: }
597: }
598: } else {
599: jac->mat = jac->pmat;
600: }
602: /* Check for null space attached to IS */
603: ilink = jac->head;
604: for (i=0; i<nsplit; i++) {
605: MatNullSpace sp;
607: PetscObjectQuery((PetscObject) ilink->is, "nullspace", (PetscObject*) &sp);
608: if (sp) {
609: MatSetNullSpace(jac->mat[i], sp);
610: }
611: ilink = ilink->next;
612: }
614: if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR) {
615: /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */
616: /* FIXME: Can/should we reuse jac->mat whenever (jac->diag_use_amat) is true? */
617: ilink = jac->head;
618: if (nsplit == 2 && jac->type == PC_COMPOSITE_MULTIPLICATIVE) {
619: /* 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 */
620: if (!jac->Afield) {
621: PetscCalloc1(nsplit,&jac->Afield);
622: MatGetSubMatrix(pc->mat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);
623: } else {
624: MatGetSubMatrix(pc->mat,ilink->next->is,ilink->is,MAT_REUSE_MATRIX,&jac->Afield[1]);
625: }
626: } else {
627: if (!jac->Afield) {
628: PetscMalloc1(nsplit,&jac->Afield);
629: for (i=0; i<nsplit; i++) {
630: MatGetSubMatrix(pc->mat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);
631: ilink = ilink->next;
632: }
633: } else {
634: for (i=0; i<nsplit; i++) {
635: MatGetSubMatrix(pc->mat,ilink->is,NULL,MAT_REUSE_MATRIX,&jac->Afield[i]);
636: ilink = ilink->next;
637: }
638: }
639: }
640: }
642: if (jac->type == PC_COMPOSITE_SCHUR) {
643: IS ccis;
644: PetscInt rstart,rend;
645: char lscname[256];
646: PetscObject LSC_L;
648: if (nsplit != 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"To use Schur complement preconditioner you must have exactly 2 fields");
650: /* When extracting off-diagonal submatrices, we take complements from this range */
651: MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);
653: /* need to handle case when one is resetting up the preconditioner */
654: if (jac->schur) {
655: KSP kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper;
657: MatSchurComplementGetKSP(jac->schur, &kspInner);
658: ilink = jac->head;
659: ISComplement(ilink->is_col,rstart,rend,&ccis);
660: if (jac->offdiag_use_amat) {
661: MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);
662: } else {
663: MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);
664: }
665: ISDestroy(&ccis);
666: ilink = ilink->next;
667: ISComplement(ilink->is_col,rstart,rend,&ccis);
668: if (jac->offdiag_use_amat) {
669: MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);
670: } else {
671: MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);
672: }
673: ISDestroy(&ccis);
674: MatSchurComplementUpdateSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);
675: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
676: MatDestroy(&jac->schurp);
677: MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);
678: }
679: if (kspA != kspInner) {
680: KSPSetOperators(kspA,jac->mat[0],jac->pmat[0]);
681: }
682: if (kspUpper != kspA) {
683: KSPSetOperators(kspUpper,jac->mat[0],jac->pmat[0]);
684: }
685: KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));
686: } else {
687: const char *Dprefix;
688: char schurprefix[256], schurmatprefix[256];
689: char schurtestoption[256];
690: MatNullSpace sp;
691: PetscBool flg;
693: /* extract the A01 and A10 matrices */
694: ilink = jac->head;
695: ISComplement(ilink->is_col,rstart,rend,&ccis);
696: if (jac->offdiag_use_amat) {
697: MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
698: } else {
699: MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
700: }
701: ISDestroy(&ccis);
702: ilink = ilink->next;
703: ISComplement(ilink->is_col,rstart,rend,&ccis);
704: if (jac->offdiag_use_amat) {
705: MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
706: } else {
707: MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
708: }
709: ISDestroy(&ccis);
711: /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */
712: MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur);
713: MatSetType(jac->schur,MATSCHURCOMPLEMENT);
714: MatSchurComplementSetSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);
715: PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
716: /* Note that the inner KSP is NOT going to inherit this prefix, and if it did, it would be reset just below. Is that what we want? */
717: MatSetOptionsPrefix(jac->schur,schurmatprefix);
718: MatSetFromOptions(jac->schur);
719: MatGetNullSpace(jac->mat[1], &sp);
720: if (sp) {
721: MatSetNullSpace(jac->schur, sp);
722: }
724: PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);
725: PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
726: if (flg) {
727: DM dmInner;
728: KSP kspInner;
730: MatSchurComplementGetKSP(jac->schur, &kspInner);
731: PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
732: /* Indent this deeper to emphasize the "inner" nature of this solver. */
733: PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject) pc, 2);
734: KSPSetOptionsPrefix(kspInner, schurprefix);
736: /* Set DM for new solver */
737: KSPGetDM(jac->head->ksp, &dmInner);
738: KSPSetDM(kspInner, dmInner);
739: KSPSetDMActive(kspInner, PETSC_FALSE);
740: } else {
741: /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or
742: * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact,
743: * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for
744: * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make
745: * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used
746: * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */
747: KSPSetType(jac->head->ksp,KSPGMRES);
748: MatSchurComplementSetKSP(jac->schur,jac->head->ksp);
749: }
750: KSPSetOperators(jac->head->ksp,jac->mat[0],jac->pmat[0]);
751: KSPSetFromOptions(jac->head->ksp);
752: MatSetFromOptions(jac->schur);
754: PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);
755: PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
756: if (flg) {
757: DM dmInner;
759: PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
760: KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper);
761: KSPSetErrorIfNotConverged(jac->kspupper,pc->erroriffailure);
762: KSPSetOptionsPrefix(jac->kspupper, schurprefix);
763: KSPGetDM(jac->head->ksp, &dmInner);
764: KSPSetDM(jac->kspupper, dmInner);
765: KSPSetDMActive(jac->kspupper, PETSC_FALSE);
766: KSPSetFromOptions(jac->kspupper);
767: KSPSetOperators(jac->kspupper,jac->mat[0],jac->pmat[0]);
768: VecDuplicate(jac->head->x, &jac->head->z);
769: } else {
770: jac->kspupper = jac->head->ksp;
771: PetscObjectReference((PetscObject) jac->head->ksp);
772: }
774: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
775: MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);
776: }
777: KSPCreate(PetscObjectComm((PetscObject)pc),&jac->kspschur);
778: KSPSetErrorIfNotConverged(jac->kspschur,pc->erroriffailure);
779: PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);
780: PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);
781: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) {
782: PC pcschur;
783: KSPGetPC(jac->kspschur,&pcschur);
784: PCSetType(pcschur,PCNONE);
785: /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */
786: } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) {
787: MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user);
788: }
789: KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));
790: KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix);
791: KSPSetOptionsPrefix(jac->kspschur, Dprefix);
792: /* propogate DM */
793: {
794: DM sdm;
795: KSPGetDM(jac->head->next->ksp, &sdm);
796: if (sdm) {
797: KSPSetDM(jac->kspschur, sdm);
798: KSPSetDMActive(jac->kspschur, PETSC_FALSE);
799: }
800: }
801: /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
802: /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */
803: KSPSetFromOptions(jac->kspschur);
804: }
806: /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */
807: PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_L",ilink->splitname);
808: PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);
809: if (!LSC_L) {PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);}
810: if (LSC_L) {PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L);}
811: PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_Lp",ilink->splitname);
812: PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);
813: if (!LSC_L) {PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);}
814: if (LSC_L) {PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L);}
815: } else {
816: /* set up the individual splits' PCs */
817: i = 0;
818: ilink = jac->head;
819: while (ilink) {
820: KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i]);
821: /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
822: if (!jac->suboptionsset) {KSPSetFromOptions(ilink->ksp);}
823: i++;
824: ilink = ilink->next;
825: }
826: }
828: jac->suboptionsset = PETSC_TRUE;
829: return(0);
830: }
832: #define FieldSplitSplitSolveAdd(ilink,xx,yy) \
833: (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
834: VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
835: KSPSolve(ilink->ksp,ilink->x,ilink->y) || \
836: VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \
837: VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE))
841: static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y)
842: {
843: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
844: PetscErrorCode ierr;
845: PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
846: KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;
849: switch (jac->schurfactorization) {
850: case PC_FIELDSPLIT_SCHUR_FACT_DIAG:
851: /* [A00 0; 0 -S], positive definite, suitable for MINRES */
852: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
853: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
854: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
855: KSPSolve(kspA,ilinkA->x,ilinkA->y);
856: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
857: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
858: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
859: VecScale(ilinkD->y,-1.);
860: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
861: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
862: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
863: break;
864: case PC_FIELDSPLIT_SCHUR_FACT_LOWER:
865: /* [A00 0; A10 S], suitable for left preconditioning */
866: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
867: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
868: KSPSolve(kspA,ilinkA->x,ilinkA->y);
869: MatMult(jac->C,ilinkA->y,ilinkD->x);
870: VecScale(ilinkD->x,-1.);
871: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
872: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
873: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
874: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
875: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
876: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
877: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
878: break;
879: case PC_FIELDSPLIT_SCHUR_FACT_UPPER:
880: /* [A00 A01; 0 S], suitable for right preconditioning */
881: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
882: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
883: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
884: MatMult(jac->B,ilinkD->y,ilinkA->x);
885: VecScale(ilinkA->x,-1.);
886: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);
887: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
888: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);
889: KSPSolve(kspA,ilinkA->x,ilinkA->y);
890: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
891: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
892: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
893: break;
894: case PC_FIELDSPLIT_SCHUR_FACT_FULL:
895: /* [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 */
896: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
897: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
898: KSPSolve(kspLower,ilinkA->x,ilinkA->y);
899: MatMult(jac->C,ilinkA->y,ilinkD->x);
900: VecScale(ilinkD->x,-1.0);
901: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
902: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
904: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
905: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
906: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
908: if (kspUpper == kspA) {
909: MatMult(jac->B,ilinkD->y,ilinkA->y);
910: VecAXPY(ilinkA->x,-1.0,ilinkA->y);
911: KSPSolve(kspA,ilinkA->x,ilinkA->y);
912: } else {
913: KSPSolve(kspA,ilinkA->x,ilinkA->y);
914: MatMult(jac->B,ilinkD->y,ilinkA->x);
915: KSPSolve(kspUpper,ilinkA->x,ilinkA->z);
916: VecAXPY(ilinkA->y,-1.0,ilinkA->z);
917: }
918: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
919: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
920: }
921: return(0);
922: }
926: static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y)
927: {
928: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
929: PetscErrorCode ierr;
930: PC_FieldSplitLink ilink = jac->head;
931: PetscInt cnt,bs;
934: if (jac->type == PC_COMPOSITE_ADDITIVE) {
935: if (jac->defaultsplit) {
936: VecGetBlockSize(x,&bs);
937: 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);
938: VecGetBlockSize(y,&bs);
939: 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);
940: VecStrideGatherAll(x,jac->x,INSERT_VALUES);
941: while (ilink) {
942: KSPSolve(ilink->ksp,ilink->x,ilink->y);
943: ilink = ilink->next;
944: }
945: VecStrideScatterAll(jac->y,y,INSERT_VALUES);
946: } else {
947: VecSet(y,0.0);
948: while (ilink) {
949: FieldSplitSplitSolveAdd(ilink,x,y);
950: ilink = ilink->next;
951: }
952: }
953: } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) {
954: VecSet(y,0.0);
955: /* solve on first block for first block variables */
956: VecScatterBegin(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);
957: VecScatterEnd(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);
958: KSPSolve(ilink->ksp,ilink->x,ilink->y);
959: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
960: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
962: /* compute the residual only onto second block variables using first block variables */
963: MatMult(jac->Afield[1],ilink->y,ilink->next->x);
964: ilink = ilink->next;
965: VecScale(ilink->x,-1.0);
966: VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
967: VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
969: /* solve on second block variables */
970: KSPSolve(ilink->ksp,ilink->x,ilink->y);
971: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
972: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
973: } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
974: if (!jac->w1) {
975: VecDuplicate(x,&jac->w1);
976: VecDuplicate(x,&jac->w2);
977: }
978: VecSet(y,0.0);
979: FieldSplitSplitSolveAdd(ilink,x,y);
980: cnt = 1;
981: while (ilink->next) {
982: ilink = ilink->next;
983: /* compute the residual only over the part of the vector needed */
984: MatMult(jac->Afield[cnt++],y,ilink->x);
985: VecScale(ilink->x,-1.0);
986: VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
987: VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
988: KSPSolve(ilink->ksp,ilink->x,ilink->y);
989: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
990: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
991: }
992: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
993: cnt -= 2;
994: while (ilink->previous) {
995: ilink = ilink->previous;
996: /* compute the residual only over the part of the vector needed */
997: MatMult(jac->Afield[cnt--],y,ilink->x);
998: VecScale(ilink->x,-1.0);
999: VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1000: VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1001: KSPSolve(ilink->ksp,ilink->x,ilink->y);
1002: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1003: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1004: }
1005: }
1006: } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type);
1007: return(0);
1008: }
1010: #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \
1011: (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
1012: VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
1013: KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \
1014: VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \
1015: VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE))
1019: static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y)
1020: {
1021: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1022: PetscErrorCode ierr;
1023: PC_FieldSplitLink ilink = jac->head;
1024: PetscInt bs;
1027: if (jac->type == PC_COMPOSITE_ADDITIVE) {
1028: if (jac->defaultsplit) {
1029: VecGetBlockSize(x,&bs);
1030: 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);
1031: VecGetBlockSize(y,&bs);
1032: 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);
1033: VecStrideGatherAll(x,jac->x,INSERT_VALUES);
1034: while (ilink) {
1035: KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);
1036: ilink = ilink->next;
1037: }
1038: VecStrideScatterAll(jac->y,y,INSERT_VALUES);
1039: } else {
1040: VecSet(y,0.0);
1041: while (ilink) {
1042: FieldSplitSplitSolveAddTranspose(ilink,x,y);
1043: ilink = ilink->next;
1044: }
1045: }
1046: } else {
1047: if (!jac->w1) {
1048: VecDuplicate(x,&jac->w1);
1049: VecDuplicate(x,&jac->w2);
1050: }
1051: VecSet(y,0.0);
1052: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1053: FieldSplitSplitSolveAddTranspose(ilink,x,y);
1054: while (ilink->next) {
1055: ilink = ilink->next;
1056: MatMultTranspose(pc->mat,y,jac->w1);
1057: VecWAXPY(jac->w2,-1.0,jac->w1,x);
1058: FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1059: }
1060: while (ilink->previous) {
1061: ilink = ilink->previous;
1062: MatMultTranspose(pc->mat,y,jac->w1);
1063: VecWAXPY(jac->w2,-1.0,jac->w1,x);
1064: FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1065: }
1066: } else {
1067: while (ilink->next) { /* get to last entry in linked list */
1068: ilink = ilink->next;
1069: }
1070: FieldSplitSplitSolveAddTranspose(ilink,x,y);
1071: while (ilink->previous) {
1072: ilink = ilink->previous;
1073: MatMultTranspose(pc->mat,y,jac->w1);
1074: VecWAXPY(jac->w2,-1.0,jac->w1,x);
1075: FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1076: }
1077: }
1078: }
1079: return(0);
1080: }
1084: static PetscErrorCode PCReset_FieldSplit(PC pc)
1085: {
1086: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1087: PetscErrorCode ierr;
1088: PC_FieldSplitLink ilink = jac->head,next;
1091: while (ilink) {
1092: KSPReset(ilink->ksp);
1093: VecDestroy(&ilink->x);
1094: VecDestroy(&ilink->y);
1095: VecDestroy(&ilink->z);
1096: VecScatterDestroy(&ilink->sctx);
1097: ISDestroy(&ilink->is);
1098: ISDestroy(&ilink->is_col);
1099: next = ilink->next;
1100: ilink = next;
1101: }
1102: PetscFree2(jac->x,jac->y);
1103: if (jac->mat && jac->mat != jac->pmat) {
1104: MatDestroyMatrices(jac->nsplits,&jac->mat);
1105: } else if (jac->mat) {
1106: jac->mat = NULL;
1107: }
1108: if (jac->pmat) {MatDestroyMatrices(jac->nsplits,&jac->pmat);}
1109: if (jac->Afield) {MatDestroyMatrices(jac->nsplits,&jac->Afield);}
1110: VecDestroy(&jac->w1);
1111: VecDestroy(&jac->w2);
1112: MatDestroy(&jac->schur);
1113: MatDestroy(&jac->schurp);
1114: MatDestroy(&jac->schur_user);
1115: KSPDestroy(&jac->kspschur);
1116: KSPDestroy(&jac->kspupper);
1117: MatDestroy(&jac->B);
1118: MatDestroy(&jac->C);
1119: jac->reset = PETSC_TRUE;
1120: return(0);
1121: }
1125: static PetscErrorCode PCDestroy_FieldSplit(PC pc)
1126: {
1127: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1128: PetscErrorCode ierr;
1129: PC_FieldSplitLink ilink = jac->head,next;
1132: PCReset_FieldSplit(pc);
1133: while (ilink) {
1134: KSPDestroy(&ilink->ksp);
1135: next = ilink->next;
1136: PetscFree(ilink->splitname);
1137: PetscFree(ilink->fields);
1138: PetscFree(ilink->fields_col);
1139: PetscFree(ilink);
1140: ilink = next;
1141: }
1142: PetscFree2(jac->x,jac->y);
1143: PetscFree(pc->data);
1144: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL);
1145: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL);
1146: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL);
1147: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL);
1148: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL);
1149: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL);
1150: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL);
1151: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL);
1152: return(0);
1153: }
1157: static PetscErrorCode PCSetFromOptions_FieldSplit(PetscOptions *PetscOptionsObject,PC pc)
1158: {
1159: PetscErrorCode ierr;
1160: PetscInt bs;
1161: PetscBool flg,stokes = PETSC_FALSE;
1162: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1163: PCCompositeType ctype;
1166: PetscOptionsHead(PetscOptionsObject,"FieldSplit options");
1167: PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL);
1168: PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);
1169: if (flg) {
1170: PCFieldSplitSetBlockSize(pc,bs);
1171: }
1172: jac->diag_use_amat = pc->useAmat;
1173: 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);
1174: jac->offdiag_use_amat = pc->useAmat;
1175: 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);
1176: /* FIXME: No programmatic equivalent to the following. */
1177: PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);
1178: if (stokes) {
1179: PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);
1180: jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF;
1181: }
1183: PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);
1184: if (flg) {
1185: PCFieldSplitSetType(pc,ctype);
1186: }
1187: /* Only setup fields once */
1188: if ((jac->bs > 0) && (jac->nsplits == 0)) {
1189: /* only allow user to set fields from command line if bs is already known.
1190: otherwise user can set them in PCFieldSplitSetDefaults() */
1191: PCFieldSplitSetRuntimeSplits_Private(pc);
1192: if (jac->splitdefined) {PetscInfo(pc,"Splits defined using the options database\n");}
1193: }
1194: if (jac->type == PC_COMPOSITE_SCHUR) {
1195: PetscOptionsGetEnum(((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);
1196: if (flg) {PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");}
1197: 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);
1198: PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSetSchurPre",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL);
1199: }
1200: PetscOptionsTail();
1201: return(0);
1202: }
1204: /*------------------------------------------------------------------------------------*/
1208: static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1209: {
1210: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1211: PetscErrorCode ierr;
1212: PC_FieldSplitLink ilink,next = jac->head;
1213: char prefix[128];
1214: PetscInt i;
1217: if (jac->splitdefined) {
1218: PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1219: return(0);
1220: }
1221: for (i=0; i<n; i++) {
1222: if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs);
1223: if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]);
1224: }
1225: PetscNew(&ilink);
1226: if (splitname) {
1227: PetscStrallocpy(splitname,&ilink->splitname);
1228: } else {
1229: PetscMalloc1(3,&ilink->splitname);
1230: PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);
1231: }
1232: PetscMalloc1(n,&ilink->fields);
1233: PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));
1234: PetscMalloc1(n,&ilink->fields_col);
1235: PetscMemcpy(ilink->fields_col,fields_col,n*sizeof(PetscInt));
1237: ilink->nfields = n;
1238: ilink->next = NULL;
1239: KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1240: KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);
1241: PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1242: KSPSetType(ilink->ksp,KSPPREONLY);
1243: PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);
1245: PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);
1246: KSPSetOptionsPrefix(ilink->ksp,prefix);
1248: if (!next) {
1249: jac->head = ilink;
1250: ilink->previous = NULL;
1251: } else {
1252: while (next->next) {
1253: next = next->next;
1254: }
1255: next->next = ilink;
1256: ilink->previous = next;
1257: }
1258: jac->nsplits++;
1259: return(0);
1260: }
1264: static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp)
1265: {
1266: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1270: PetscMalloc1(jac->nsplits,subksp);
1271: MatSchurComplementGetKSP(jac->schur,*subksp);
1273: (*subksp)[1] = jac->kspschur;
1274: if (n) *n = jac->nsplits;
1275: return(0);
1276: }
1280: static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp)
1281: {
1282: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1283: PetscErrorCode ierr;
1284: PetscInt cnt = 0;
1285: PC_FieldSplitLink ilink = jac->head;
1288: PetscMalloc1(jac->nsplits,subksp);
1289: while (ilink) {
1290: (*subksp)[cnt++] = ilink->ksp;
1291: ilink = ilink->next;
1292: }
1293: 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);
1294: if (n) *n = jac->nsplits;
1295: return(0);
1296: }
1300: static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is)
1301: {
1302: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1303: PetscErrorCode ierr;
1304: PC_FieldSplitLink ilink, next = jac->head;
1305: char prefix[128];
1308: if (jac->splitdefined) {
1309: PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1310: return(0);
1311: }
1312: PetscNew(&ilink);
1313: if (splitname) {
1314: PetscStrallocpy(splitname,&ilink->splitname);
1315: } else {
1316: PetscMalloc1(8,&ilink->splitname);
1317: PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);
1318: }
1319: PetscObjectReference((PetscObject)is);
1320: ISDestroy(&ilink->is);
1321: ilink->is = is;
1322: PetscObjectReference((PetscObject)is);
1323: ISDestroy(&ilink->is_col);
1324: ilink->is_col = is;
1325: ilink->next = NULL;
1326: KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1327: KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);
1328: PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1329: KSPSetType(ilink->ksp,KSPPREONLY);
1330: PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);
1332: PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);
1333: KSPSetOptionsPrefix(ilink->ksp,prefix);
1335: if (!next) {
1336: jac->head = ilink;
1337: ilink->previous = NULL;
1338: } else {
1339: while (next->next) {
1340: next = next->next;
1341: }
1342: next->next = ilink;
1343: ilink->previous = next;
1344: }
1345: jac->nsplits++;
1346: return(0);
1347: }
1351: /*@
1352: PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner
1354: Logically Collective on PC
1356: Input Parameters:
1357: + pc - the preconditioner context
1358: . splitname - name of this split, if NULL the number of the split is used
1359: . n - the number of fields in this split
1360: - fields - the fields in this split
1362: Level: intermediate
1364: Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field.
1366: The PCFieldSplitSetFields() is for defining fields as strided blocks. For example, if the block
1367: 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
1368: 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x....
1369: where the numbered entries indicate what is in the field.
1371: This function is called once per split (it creates a new split each time). Solve options
1372: for this split will be available under the prefix -fieldsplit_SPLITNAME_.
1374: Developer Note: This routine does not actually create the IS representing the split, that is delayed
1375: until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be
1376: available when this routine is called.
1378: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS()
1380: @*/
1381: PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1382: {
1388: if (n < 1) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname);
1390: PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt*,const PetscInt*),(pc,splitname,n,fields,fields_col));
1391: return(0);
1392: }
1396: /*@
1397: PCFieldSplitSetDiagUseAmat - set flag indicating whether to extract diagonal blocks from Amat (rather than Pmat)
1399: Logically Collective on PC
1401: Input Parameters:
1402: + pc - the preconditioner object
1403: - flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from
1405: Options Database:
1406: . -pc_fieldsplit_diag_use_amat
1408: Level: intermediate
1410: .seealso: PCFieldSplitGetDiagUseAmat(), PCFieldSplitSetOffDiagUseAmat(), PCFIELDSPLIT
1412: @*/
1413: PetscErrorCode PCFieldSplitSetDiagUseAmat(PC pc,PetscBool flg)
1414: {
1415: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1416: PetscBool isfs;
1421: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1422: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1423: jac->diag_use_amat = flg;
1424: return(0);
1425: }
1429: /*@
1430: PCFieldSplitGetDiagUseAmat - get the flag indicating whether to extract diagonal blocks from Amat (rather than Pmat)
1432: Logically Collective on PC
1434: Input Parameters:
1435: . pc - the preconditioner object
1437: Output Parameters:
1438: . flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from
1441: Level: intermediate
1443: .seealso: PCFieldSplitSetDiagUseAmat(), PCFieldSplitGetOffDiagUseAmat(), PCFIELDSPLIT
1445: @*/
1446: PetscErrorCode PCFieldSplitGetDiagUseAmat(PC pc,PetscBool *flg)
1447: {
1448: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1449: PetscBool isfs;
1455: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1456: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1457: *flg = jac->diag_use_amat;
1458: return(0);
1459: }
1463: /*@
1464: PCFieldSplitSetOffDiagUseAmat - set flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat)
1466: Logically Collective on PC
1468: Input Parameters:
1469: + pc - the preconditioner object
1470: - flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from
1472: Options Database:
1473: . -pc_fieldsplit_off_diag_use_amat
1475: Level: intermediate
1477: .seealso: PCFieldSplitGetOffDiagUseAmat(), PCFieldSplitSetDiagUseAmat(), PCFIELDSPLIT
1479: @*/
1480: PetscErrorCode PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg)
1481: {
1482: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1483: PetscBool isfs;
1488: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1489: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1490: jac->offdiag_use_amat = flg;
1491: return(0);
1492: }
1496: /*@
1497: PCFieldSplitGetOffDiagUseAmat - get the flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat)
1499: Logically Collective on PC
1501: Input Parameters:
1502: . pc - the preconditioner object
1504: Output Parameters:
1505: . flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from
1508: Level: intermediate
1510: .seealso: PCFieldSplitSetOffDiagUseAmat(), PCFieldSplitGetDiagUseAmat(), PCFIELDSPLIT
1512: @*/
1513: PetscErrorCode PCFieldSplitGetOffDiagUseAmat(PC pc,PetscBool *flg)
1514: {
1515: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1516: PetscBool isfs;
1522: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1523: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1524: *flg = jac->offdiag_use_amat;
1525: return(0);
1526: }
1532: /*@C
1533: PCFieldSplitSetIS - Sets the exact elements for field
1535: Logically Collective on PC
1537: Input Parameters:
1538: + pc - the preconditioner context
1539: . splitname - name of this split, if NULL the number of the split is used
1540: - is - the index set that defines the vector elements in this field
1543: Notes:
1544: Use PCFieldSplitSetFields(), for fields defined by strided types.
1546: This function is called once per split (it creates a new split each time). Solve options
1547: for this split will be available under the prefix -fieldsplit_SPLITNAME_.
1549: Level: intermediate
1551: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize()
1553: @*/
1554: PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is)
1555: {
1562: PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));
1563: return(0);
1564: }
1568: /*@
1569: PCFieldSplitGetIS - Retrieves the elements for a field as an IS
1571: Logically Collective on PC
1573: Input Parameters:
1574: + pc - the preconditioner context
1575: - splitname - name of this split
1577: Output Parameter:
1578: - is - the index set that defines the vector elements in this field, or NULL if the field is not found
1580: Level: intermediate
1582: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS()
1584: @*/
1585: PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is)
1586: {
1593: {
1594: PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;
1595: PC_FieldSplitLink ilink = jac->head;
1596: PetscBool found;
1598: *is = NULL;
1599: while (ilink) {
1600: PetscStrcmp(ilink->splitname, splitname, &found);
1601: if (found) {
1602: *is = ilink->is;
1603: break;
1604: }
1605: ilink = ilink->next;
1606: }
1607: }
1608: return(0);
1609: }
1613: /*@
1614: PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the
1615: fieldsplit preconditioner. If not set the matrix block size is used.
1617: Logically Collective on PC
1619: Input Parameters:
1620: + pc - the preconditioner context
1621: - bs - the block size
1623: Level: intermediate
1625: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields()
1627: @*/
1628: PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs)
1629: {
1635: PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));
1636: return(0);
1637: }
1641: /*@C
1642: PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits
1644: Collective on KSP
1646: Input Parameter:
1647: . pc - the preconditioner context
1649: Output Parameters:
1650: + n - the number of splits
1651: - pc - the array of KSP contexts
1653: Note:
1654: After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user
1655: (not the KSP just the array that contains them).
1657: You must call KSPSetUp() before calling PCFieldSplitGetSubKSP().
1659: Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs.
1660: You can call PCFieldSplitGetSubKSP(pc,n,NULL_OBJECT,ierr) to determine how large the
1661: KSP array must be.
1664: Level: advanced
1666: .seealso: PCFIELDSPLIT
1667: @*/
1668: PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
1669: {
1675: PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));
1676: return(0);
1677: }
1681: /*@
1682: PCFieldSplitSetSchurPre - Indicates if the Schur complement is preconditioned by a preconditioner constructed by the
1683: A11 matrix. Otherwise no preconditioner is used.
1685: Collective on PC
1687: Input Parameters:
1688: + pc - the preconditioner context
1689: . ptype - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11 (default), PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_PRE_USER
1690: - userpre - matrix to use for preconditioning, or NULL
1692: Options Database:
1693: . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> default is a11
1695: Notes:
1696: $ If ptype is
1697: $ user then the preconditioner for the Schur complement is generated by the provided matrix (pre argument
1698: $ to this function).
1699: $ a11 then the preconditioner for the Schur complement is generated by the block diagonal part of the preconditioner
1700: $ matrix associated with the Schur complement (i.e. A11), not he Schur complement matrix
1701: $ full then the preconditioner uses the exact Schur complement (this is expensive)
1702: $ self the preconditioner for the Schur complement is generated from the Schur complement matrix itself:
1703: $ The only preconditioner that currently works directly with the Schur complement matrix object is the PCLSC
1704: $ preconditioner
1705: $ selfp then the preconditioning matrix is an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01
1706: $ This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be
1707: $ lumped before extracting the diagonal: -fieldsplit_1_mat_schur_complement_ainv_type lump; diag is the default.
1709: When solving a saddle point problem, where the A11 block is identically zero, using a11 as the ptype only makes sense
1710: with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and
1711: -fieldsplit_1_pc_type lsc which uses the least squares commutator to compute a preconditioner for the Schur complement.
1713: Level: intermediate
1715: .seealso: PCFieldSplitGetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType,
1716: MatSchurComplementSetAinvType(), PCLSC
1718: @*/
1719: PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
1720: {
1725: PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));
1726: return(0);
1727: }
1728: PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) {return PCFieldSplitSetSchurPre(pc,ptype,pre);} /* Deprecated name */
1732: /*@
1733: PCFieldSplitGetSchurPre - For Schur complement fieldsplit, determine how the Schur complement will be
1734: preconditioned. See PCFieldSplitSetSchurPre() for details.
1736: Logically Collective on PC
1738: Input Parameters:
1739: . pc - the preconditioner context
1741: Output Parameters:
1742: + ptype - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11, PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_PRE_USER
1743: - userpre - matrix to use for preconditioning (with PC_FIELDSPLIT_PRE_USER), or NULL
1745: Level: intermediate
1747: .seealso: PCFieldSplitSetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC
1749: @*/
1750: PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
1751: {
1756: PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre));
1757: return(0);
1758: }
1762: /*@
1763: PCFieldSplitSchurGetS - extract the MatSchurComplement object used by this PC in case it needs to be configured separately
1765: Not collective
1767: Input Parameter:
1768: . pc - the preconditioner context
1770: Output Parameter:
1771: . S - the Schur complement matrix
1773: Notes:
1774: This matrix should not be destroyed using MatDestroy(); rather, use PCFieldSplitSchurRestoreS().
1776: Level: advanced
1778: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurRestoreS()
1780: @*/
1781: PetscErrorCode PCFieldSplitSchurGetS(PC pc,Mat *S)
1782: {
1784: const char* t;
1785: PetscBool isfs;
1786: PC_FieldSplit *jac;
1790: PetscObjectGetType((PetscObject)pc,&t);
1791: PetscStrcmp(t,PCFIELDSPLIT,&isfs);
1792: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
1793: jac = (PC_FieldSplit*)pc->data;
1794: if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
1795: if (S) *S = jac->schur;
1796: return(0);
1797: }
1801: /*@
1802: PCFieldSplitSchurRestoreS - restores the MatSchurComplement object used by this PC
1804: Not collective
1806: Input Parameters:
1807: + pc - the preconditioner context
1808: . S - the Schur complement matrix
1810: Level: advanced
1812: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurGetS()
1814: @*/
1815: PetscErrorCode PCFieldSplitSchurRestoreS(PC pc,Mat *S)
1816: {
1818: const char* t;
1819: PetscBool isfs;
1820: PC_FieldSplit *jac;
1824: PetscObjectGetType((PetscObject)pc,&t);
1825: PetscStrcmp(t,PCFIELDSPLIT,&isfs);
1826: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
1827: jac = (PC_FieldSplit*)pc->data;
1828: if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
1829: if (!S || *S != jac->schur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten");
1830: return(0);
1831: }
1836: static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
1837: {
1838: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1842: jac->schurpre = ptype;
1843: if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) {
1844: MatDestroy(&jac->schur_user);
1845: jac->schur_user = pre;
1846: PetscObjectReference((PetscObject)jac->schur_user);
1847: }
1848: return(0);
1849: }
1853: static PetscErrorCode PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
1854: {
1855: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1858: *ptype = jac->schurpre;
1859: *pre = jac->schur_user;
1860: return(0);
1861: }
1865: /*@
1866: PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain
1868: Collective on PC
1870: Input Parameters:
1871: + pc - the preconditioner context
1872: - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default
1874: Options Database:
1875: . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full
1878: Level: intermediate
1880: Notes:
1881: The FULL factorization is
1883: $ (A B) = (1 0) (A 0) (1 Ainv*B)
1884: $ (C D) (C*Ainv 1) (0 S) (0 1 )
1886: where S = D - C*Ainv*B. In practice, the full factorization is applied via block triangular solves with the grouping L*(D*U). UPPER uses D*U, LOWER uses L*D,
1887: 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).
1889: If applied exactly, FULL factorization is a direct solver. The preconditioned operator with LOWER or UPPER has all eigenvalues equal to 1 and minimal polynomial
1890: of degree 2, so KSPGMRES converges in 2 iterations. If the iteration count is very low, consider using KSPFGMRES or KSPGCR which can use one less preconditioner
1891: application in this case. Note that the preconditioned operator may be highly non-normal, so such fast convergence may not be observed in practice. With DIAG,
1892: the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so KSPGMRES converges in at most 4 iterations.
1894: For symmetric problems in which A is positive definite and S is negative definite, DIAG can be used with KSPMINRES. Note that a flexible method like KSPFGMRES
1895: or KSPGCR must be used if the fieldsplit preconditioner is nonlinear (e.g. a few iterations of a Krylov method is used inside a split).
1897: References:
1898: Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000) pp. 1969-1972.
1900: Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001), pp. 1050-1051.
1902: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType
1903: @*/
1904: PetscErrorCode PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype)
1905: {
1910: PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));
1911: return(0);
1912: }
1916: static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype)
1917: {
1918: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1921: jac->schurfactorization = ftype;
1922: return(0);
1923: }
1927: /*@C
1928: PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement
1930: Collective on KSP
1932: Input Parameter:
1933: . pc - the preconditioner context
1935: Output Parameters:
1936: + A00 - the (0,0) block
1937: . A01 - the (0,1) block
1938: . A10 - the (1,0) block
1939: - A11 - the (1,1) block
1941: Level: advanced
1943: .seealso: PCFIELDSPLIT
1944: @*/
1945: PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11)
1946: {
1947: PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;
1951: if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach.");
1952: if (A00) *A00 = jac->pmat[0];
1953: if (A01) *A01 = jac->B;
1954: if (A10) *A10 = jac->C;
1955: if (A11) *A11 = jac->pmat[1];
1956: return(0);
1957: }
1961: static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type)
1962: {
1963: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1967: jac->type = type;
1968: if (type == PC_COMPOSITE_SCHUR) {
1969: pc->ops->apply = PCApply_FieldSplit_Schur;
1970: pc->ops->view = PCView_FieldSplit_Schur;
1972: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);
1973: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit);
1974: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit);
1975: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);
1977: } else {
1978: pc->ops->apply = PCApply_FieldSplit;
1979: pc->ops->view = PCView_FieldSplit;
1981: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
1982: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",0);
1983: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",0);
1984: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);
1985: }
1986: return(0);
1987: }
1991: static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs)
1992: {
1993: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1996: if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs);
1997: 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);
1998: jac->bs = bs;
1999: return(0);
2000: }
2004: /*@
2005: PCFieldSplitSetType - Sets the type of fieldsplit preconditioner.
2007: Collective on PC
2009: Input Parameter:
2010: . pc - the preconditioner context
2011: . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR
2013: Options Database Key:
2014: . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type
2016: Level: Intermediate
2018: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
2020: .seealso: PCCompositeSetType()
2022: @*/
2023: PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type)
2024: {
2029: PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));
2030: return(0);
2031: }
2035: /*@
2036: PCFieldSplitGetType - Gets the type of fieldsplit preconditioner.
2038: Not collective
2040: Input Parameter:
2041: . pc - the preconditioner context
2043: Output Parameter:
2044: . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR
2046: Level: Intermediate
2048: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
2049: .seealso: PCCompositeSetType()
2050: @*/
2051: PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type)
2052: {
2053: PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;
2058: *type = jac->type;
2059: return(0);
2060: }
2064: /*@
2065: PCFieldSplitSetDMSplits - Flags whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible.
2067: Logically Collective
2069: Input Parameters:
2070: + pc - the preconditioner context
2071: - flg - boolean indicating whether to use field splits defined by the DM
2073: Options Database Key:
2074: . -pc_fieldsplit_dm_splits
2076: Level: Intermediate
2078: .keywords: PC, DM, composite preconditioner, additive, multiplicative
2080: .seealso: PCFieldSplitGetDMSplits()
2082: @*/
2083: PetscErrorCode PCFieldSplitSetDMSplits(PC pc,PetscBool flg)
2084: {
2085: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2086: PetscBool isfs;
2092: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2093: if (isfs) {
2094: jac->dm_splits = flg;
2095: }
2096: return(0);
2097: }
2102: /*@
2103: PCFieldSplitGetDMSplits - Returns flag indicating whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible.
2105: Logically Collective
2107: Input Parameter:
2108: . pc - the preconditioner context
2110: Output Parameter:
2111: . flg - boolean indicating whether to use field splits defined by the DM
2113: Level: Intermediate
2115: .keywords: PC, DM, composite preconditioner, additive, multiplicative
2117: .seealso: PCFieldSplitSetDMSplits()
2119: @*/
2120: PetscErrorCode PCFieldSplitGetDMSplits(PC pc,PetscBool* flg)
2121: {
2122: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2123: PetscBool isfs;
2129: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2130: if (isfs) {
2131: if(flg) *flg = jac->dm_splits;
2132: }
2133: return(0);
2134: }
2136: /* -------------------------------------------------------------------------------------*/
2137: /*MC
2138: PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual
2139: fields or groups of fields. See the users manual section "Solving Block Matrices" for more details.
2141: To set options on the solvers for each block append -fieldsplit_ to all the PC
2142: options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1
2144: To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP()
2145: and set the options directly on the resulting KSP object
2147: Level: intermediate
2149: Options Database Keys:
2150: + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split
2151: . -pc_fieldsplit_default - automatically add any fields to additional splits that have not
2152: been supplied explicitly by -pc_fieldsplit_%d_fields
2153: . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields)
2154: . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur> - type of relaxation or factorization splitting
2155: . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11
2156: . -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero or negative diagonal and uses Schur complement with no preconditioner as the solver
2158: - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_
2159: for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields
2161: Notes:
2162: Use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS()
2163: to define a field by an arbitrary collection of entries.
2165: If no fields are set the default is used. The fields are defined by entries strided by bs,
2166: beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(),
2167: if this is not called the block size defaults to the blocksize of the second matrix passed
2168: to KSPSetOperators()/PCSetOperators().
2170: $ For the Schur complement preconditioner if J = ( A00 A01 )
2171: $ ( A10 A11 )
2172: $ the preconditioner using full factorization is
2173: $ ( I -ksp(A00) A01 ) ( inv(A00) 0 ) ( I 0 )
2174: $ ( 0 I ) ( 0 ksp(S) ) ( -A10 ksp(A00) I )
2175: where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. S is the Schur complement
2176: $ S = A11 - A10 ksp(A00) A01
2177: 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
2178: in providing the SECOND split or 1 if not give). For PCFieldSplitGetKSP() when field number is 0,
2179: it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default
2180: A11 is used to construct a preconditioner for S, use PCFieldSplitSetSchurPre() to turn on or off this
2181: option. You can use the preconditioner PCLSC to precondition the Schur complement with -fieldsplit_1_pc_type lsc.
2182: When option -fieldsplit_schur_precondition selfp is given, an approximation to S is assembled --
2183: Sp = A11 - A10 inv(diag(A00)) A01, which has type AIJ and can be used with a variety of preconditioners
2184: (e.g., -fieldsplit_1_pc_type asm). Optionally, A00 can be lumped before extracting the diagonal:
2185: -fieldsplit_1_mat_schur_complement_ainv_type lump; diag is the default.
2187: The factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above,
2188: diag gives
2189: $ ( inv(A00) 0 )
2190: $ ( 0 -ksp(S) )
2191: note that slightly counter intuitively there is a negative in front of the ksp(S) so that the preconditioner is positive definite. The lower factorization is the inverse of
2192: $ ( A00 0 )
2193: $ ( A10 S )
2194: where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of
2195: $ ( A00 A01 )
2196: $ ( 0 S )
2197: where again the inverses of A00 and S are applied using KSPs.
2199: If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS
2200: is used automatically for a second block.
2202: The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1.
2203: Generally it should be used with the AIJ format.
2205: The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see,
2206: for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT
2207: inside a smoother resulting in "Distributive Smoothers".
2209: Concepts: physics based preconditioners, block preconditioners
2211: There is a nice discussion of block preconditioners in
2213: [El08] A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier-Stokes equations
2214: Howard Elman, V.E. Howle, John Shadid, Robert Shuttleworth, Ray Tuminaro, Journal of Computational Physics 227 (2008) 1790--1808
2215: http://chess.cs.umd.edu/~elman/papers/tax.pdf
2217: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC,
2218: PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre(),
2219: MatSchurComplementSetAinvType()
2220: M*/
2224: PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc)
2225: {
2227: PC_FieldSplit *jac;
2230: PetscNewLog(pc,&jac);
2232: jac->bs = -1;
2233: jac->nsplits = 0;
2234: jac->type = PC_COMPOSITE_MULTIPLICATIVE;
2235: jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */
2236: jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL;
2237: jac->dm_splits = PETSC_TRUE;
2239: pc->data = (void*)jac;
2241: pc->ops->apply = PCApply_FieldSplit;
2242: pc->ops->applytranspose = PCApplyTranspose_FieldSplit;
2243: pc->ops->setup = PCSetUp_FieldSplit;
2244: pc->ops->reset = PCReset_FieldSplit;
2245: pc->ops->destroy = PCDestroy_FieldSplit;
2246: pc->ops->setfromoptions = PCSetFromOptions_FieldSplit;
2247: pc->ops->view = PCView_FieldSplit;
2248: pc->ops->applyrichardson = 0;
2250: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2251: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);
2252: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);
2253: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);
2254: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);
2255: return(0);
2256: }