Actual source code: fieldsplit.c
petsc-3.5.4 2015-05-23
2: #include <petsc-private/pcimpl.h> /*I "petscpc.h" I*/
3: #include <petscdm.h>
5: /*
6: There is a nice discussion of block preconditioners in
8: [El08] A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier-Stokes equations
9: Howard Elman, V.E. Howle, John Shadid, Robert Shuttleworth, Ray Tuminaro, Journal of Computational Physics 227 (2008) 1790--1808
10: http://chess.cs.umd.edu/~elman/papers/tax.pdf
11: */
13: const char *const PCFieldSplitSchurPreTypes[] = {"SELF","SELFP","A11","USER","FULL","PCFieldSplitSchurPreType","PC_FIELDSPLIT_SCHUR_PRE_",0};
14: const char *const PCFieldSplitSchurFactTypes[] = {"DIAG","LOWER","UPPER","FULL","PCFieldSplitSchurFactType","PC_FIELDSPLIT_SCHUR_FACT_",0};
16: typedef struct _PC_FieldSplitLink *PC_FieldSplitLink;
17: struct _PC_FieldSplitLink {
18: KSP ksp;
19: Vec x,y,z;
20: char *splitname;
21: PetscInt nfields;
22: PetscInt *fields,*fields_col;
23: VecScatter sctx;
24: IS is,is_col;
25: PC_FieldSplitLink next,previous;
26: };
28: typedef struct {
29: PCCompositeType type;
30: PetscBool defaultsplit; /* Flag for a system with a set of 'k' scalar fields with the same layout (and bs = k) */
31: PetscBool splitdefined; /* Flag is set after the splits have been defined, to prevent more splits from being added */
32: PetscInt bs; /* Block size for IS and Mat structures */
33: PetscInt nsplits; /* Number of field divisions defined */
34: Vec *x,*y,w1,w2;
35: Mat *mat; /* The diagonal block for each split */
36: Mat *pmat; /* The preconditioning diagonal block for each split */
37: Mat *Afield; /* The rows of the matrix associated with each split */
38: PetscBool issetup;
40: /* Only used when Schur complement preconditioning is used */
41: Mat B; /* The (0,1) block */
42: Mat C; /* The (1,0) block */
43: Mat schur; /* The Schur complement S = A11 - A10 A00^{-1} A01, the KSP here, kspinner, is H_1 in [El08] */
44: Mat schurp; /* Assembled approximation to S built by MatSchurComplement to be used as a preconditioning matrix when solving with S */
45: Mat schur_user; /* User-provided preconditioning matrix for the Schur complement */
46: PCFieldSplitSchurPreType schurpre; /* Determines which preconditioning matrix is used for the Schur complement */
47: PCFieldSplitSchurFactType schurfactorization;
48: KSP kspschur; /* The solver for S */
49: KSP kspupper; /* The solver for A in the upper diagonal part of the factorization (H_2 in [El08]) */
50: PC_FieldSplitLink head;
51: PetscBool reset; /* indicates PCReset() has been last called on this object, hack */
52: PetscBool suboptionsset; /* Indicates that the KSPSetFromOptions() has been called on the sub-KSPs */
53: PetscBool dm_splits; /* Whether to use DMCreateFieldDecomposition() whenever possible */
54: PetscBool diag_use_amat; /* Whether to extract diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */
55: PetscBool offdiag_use_amat; /* Whether to extract off-diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */
56: } PC_FieldSplit;
58: /*
59: Notes: there is no particular reason that pmat, x, and y are stored as arrays in PC_FieldSplit instead of
60: inside PC_FieldSplitLink, just historical. If you want to be able to add new fields after already using the
61: PC you could change this.
62: */
64: /* This helper is so that setting a user-provided preconditioning matrix is orthogonal to choosing to use it. This way the
65: * application-provided FormJacobian can provide this matrix without interfering with the user's (command-line) choices. */
66: static Mat FieldSplitSchurPre(PC_FieldSplit *jac)
67: {
68: switch (jac->schurpre) {
69: case PC_FIELDSPLIT_SCHUR_PRE_SELF: return jac->schur;
70: case PC_FIELDSPLIT_SCHUR_PRE_SELFP: return jac->schurp;
71: case PC_FIELDSPLIT_SCHUR_PRE_A11: return jac->pmat[1];
72: case PC_FIELDSPLIT_SCHUR_PRE_FULL: /* We calculate this and store it in schur_user */
73: case PC_FIELDSPLIT_SCHUR_PRE_USER: /* Use a user-provided matrix if it is given, otherwise diagonal block */
74: default:
75: return jac->schur_user ? jac->schur_user : jac->pmat[1];
76: }
77: }
80: #include <petscdraw.h>
83: static PetscErrorCode PCView_FieldSplit(PC pc,PetscViewer viewer)
84: {
85: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
86: PetscErrorCode ierr;
87: PetscBool iascii,isdraw;
88: PetscInt i,j;
89: PC_FieldSplitLink ilink = jac->head;
92: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
93: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
94: if (iascii) {
95: if (jac->bs > 0) {
96: PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs);
97: } else {
98: PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D\n",PCCompositeTypes[jac->type],jac->nsplits);
99: }
100: if (pc->useAmat) {
101: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");
102: }
103: if (jac->diag_use_amat) {
104: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for diagonal blocks\n");
105: }
106: if (jac->offdiag_use_amat) {
107: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for off-diagonal blocks\n");
108: }
109: PetscViewerASCIIPrintf(viewer," Solver info for each split is in the following KSP objects:\n");
110: PetscViewerASCIIPushTab(viewer);
111: for (i=0; i<jac->nsplits; i++) {
112: if (ilink->fields) {
113: PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);
114: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
115: for (j=0; j<ilink->nfields; j++) {
116: if (j > 0) {
117: PetscViewerASCIIPrintf(viewer,",");
118: }
119: PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);
120: }
121: PetscViewerASCIIPrintf(viewer,"\n");
122: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
123: } else {
124: PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);
125: }
126: KSPView(ilink->ksp,viewer);
127: ilink = ilink->next;
128: }
129: PetscViewerASCIIPopTab(viewer);
130: }
132: if (isdraw) {
133: PetscDraw draw;
134: PetscReal x,y,w,wd;
136: PetscViewerDrawGetDraw(viewer,0,&draw);
137: PetscDrawGetCurrentPoint(draw,&x,&y);
138: w = 2*PetscMin(1.0 - x,x);
139: wd = w/(jac->nsplits + 1);
140: x = x - wd*(jac->nsplits-1)/2.0;
141: for (i=0; i<jac->nsplits; i++) {
142: PetscDrawPushCurrentPoint(draw,x,y);
143: KSPView(ilink->ksp,viewer);
144: PetscDrawPopCurrentPoint(draw);
145: x += wd;
146: ilink = ilink->next;
147: }
148: }
149: return(0);
150: }
154: static PetscErrorCode PCView_FieldSplit_Schur(PC pc,PetscViewer viewer)
155: {
156: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
157: PetscErrorCode ierr;
158: PetscBool iascii,isdraw;
159: PetscInt i,j;
160: PC_FieldSplitLink ilink = jac->head;
163: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
164: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
165: if (iascii) {
166: if (jac->bs > 0) {
167: PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, blocksize = %D, factorization %s\n",jac->bs,PCFieldSplitSchurFactTypes[jac->schurfactorization]);
168: } else {
169: PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, factorization %s\n",PCFieldSplitSchurFactTypes[jac->schurfactorization]);
170: }
171: if (pc->useAmat) {
172: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");
173: }
174: switch (jac->schurpre) {
175: case PC_FIELDSPLIT_SCHUR_PRE_SELF:
176: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from S itself\n");break;
177: case PC_FIELDSPLIT_SCHUR_PRE_SELFP:
178: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from Sp, an assembled approximation to S, which uses (the lumped) A00's diagonal's inverse\n");break;
179: case PC_FIELDSPLIT_SCHUR_PRE_A11:
180: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n");break;
181: case PC_FIELDSPLIT_SCHUR_PRE_FULL:
182: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from the exact Schur complement\n");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: PetscDrawBoxedString(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: PetscDrawBoxedString(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: }
277: /* Precondition: jac->bs is set to a meaningful value */
278: static PetscErrorCode PCFieldSplitSetRuntimeSplits_Private(PC pc)
279: {
281: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
282: PetscInt i,nfields,*ifields,nfields_col,*ifields_col;
283: PetscBool flg,flg_col;
284: char optionname[128],splitname[8],optionname_col[128];
287: PetscMalloc1(jac->bs,&ifields);
288: PetscMalloc1(jac->bs,&ifields_col);
289: for (i=0,flg=PETSC_TRUE;; i++) {
290: PetscSNPrintf(splitname,sizeof(splitname),"%D",i);
291: PetscSNPrintf(optionname,sizeof(optionname),"-pc_fieldsplit_%D_fields",i);
292: PetscSNPrintf(optionname_col,sizeof(optionname_col),"-pc_fieldsplit_%D_fields_col",i);
293: nfields = jac->bs;
294: nfields_col = jac->bs;
295: PetscOptionsGetIntArray(((PetscObject)pc)->prefix,optionname,ifields,&nfields,&flg);
296: PetscOptionsGetIntArray(((PetscObject)pc)->prefix,optionname_col,ifields_col,&nfields_col,&flg_col);
297: if (!flg) break;
298: else if (flg && !flg_col) {
299: if (!nfields) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields");
300: PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields);
301: } else {
302: if (!nfields || !nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields");
303: if (nfields != nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Number of row and column fields must match");
304: PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields_col);
305: }
306: }
307: if (i > 0) {
308: /* Makes command-line setting of splits take precedence over setting them in code.
309: Otherwise subsequent calls to PCFieldSplitSetIS() or PCFieldSplitSetFields() would
310: create new splits, which would probably not be what the user wanted. */
311: jac->splitdefined = PETSC_TRUE;
312: }
313: PetscFree(ifields);
314: PetscFree(ifields_col);
315: return(0);
316: }
320: static PetscErrorCode PCFieldSplitSetDefaults(PC pc)
321: {
322: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
323: PetscErrorCode ierr;
324: PC_FieldSplitLink ilink = jac->head;
325: PetscBool fieldsplit_default = PETSC_FALSE,stokes = PETSC_FALSE,coupling = PETSC_FALSE;
326: PetscInt i;
329: /*
330: Kinda messy, but at least this now uses DMCreateFieldDecomposition() even with jac->reset.
331: Should probably be rewritten.
332: */
333: if (!ilink || jac->reset) {
334: PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);
335: PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_coupling",&coupling,NULL);
336: if (pc->dm && jac->dm_splits && !stokes && !coupling) {
337: PetscInt numFields, f, i, j;
338: char **fieldNames;
339: IS *fields;
340: DM *dms;
341: DM subdm[128];
342: PetscBool flg;
344: DMCreateFieldDecomposition(pc->dm, &numFields, &fieldNames, &fields, &dms);
345: /* Allow the user to prescribe the splits */
346: for (i = 0, flg = PETSC_TRUE;; i++) {
347: PetscInt ifields[128];
348: IS compField;
349: char optionname[128], splitname[8];
350: PetscInt nfields = numFields;
352: PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%D_fields", i);
353: PetscOptionsGetIntArray(((PetscObject) pc)->prefix, optionname, ifields, &nfields, &flg);
354: if (!flg) break;
355: if (numFields > 128) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot currently support %d > 128 fields", numFields);
356: DMCreateSubDM(pc->dm, nfields, ifields, &compField, &subdm[i]);
357: if (nfields == 1) {
358: PCFieldSplitSetIS(pc, fieldNames[ifields[0]], compField);
359: /* PetscPrintf(PetscObjectComm((PetscObject)pc), "%s Field Indices:", fieldNames[ifields[0]]);
360: ISView(compField, NULL); */
361: } else {
362: PetscSNPrintf(splitname, sizeof(splitname), "%D", i);
363: PCFieldSplitSetIS(pc, splitname, compField);
364: /* PetscPrintf(PetscObjectComm((PetscObject)pc), "%s Field Indices:", splitname);
365: ISView(compField, NULL); */
366: }
367: ISDestroy(&compField);
368: for (j = 0; j < nfields; ++j) {
369: f = ifields[j];
370: PetscFree(fieldNames[f]);
371: ISDestroy(&fields[f]);
372: }
373: }
374: if (i == 0) {
375: for (f = 0; f < numFields; ++f) {
376: PCFieldSplitSetIS(pc, fieldNames[f], fields[f]);
377: PetscFree(fieldNames[f]);
378: ISDestroy(&fields[f]);
379: }
380: } else {
381: for (j=0; j<numFields; j++) {
382: DMDestroy(dms+j);
383: }
384: PetscFree(dms);
385: PetscMalloc1(i, &dms);
386: for (j = 0; j < i; ++j) dms[j] = subdm[j];
387: }
388: PetscFree(fieldNames);
389: PetscFree(fields);
390: if (dms) {
391: PetscInfo(pc, "Setting up physics based fieldsplit preconditioner using the embedded DM\n");
392: for (ilink = jac->head, i = 0; ilink; ilink = ilink->next, ++i) {
393: const char *prefix;
394: PetscObjectGetOptionsPrefix((PetscObject)(ilink->ksp),&prefix);
395: PetscObjectSetOptionsPrefix((PetscObject)(dms[i]), prefix);
396: KSPSetDM(ilink->ksp, dms[i]);
397: KSPSetDMActive(ilink->ksp, PETSC_FALSE);
398: PetscObjectIncrementTabLevel((PetscObject)dms[i],(PetscObject)ilink->ksp,0);
399: DMDestroy(&dms[i]);
400: }
401: PetscFree(dms);
402: }
403: } else {
404: if (jac->bs <= 0) {
405: if (pc->pmat) {
406: MatGetBlockSize(pc->pmat,&jac->bs);
407: } else jac->bs = 1;
408: }
410: if (stokes) {
411: IS zerodiags,rest;
412: PetscInt nmin,nmax;
414: MatGetOwnershipRange(pc->mat,&nmin,&nmax);
415: MatFindZeroDiagonals(pc->mat,&zerodiags);
416: ISComplement(zerodiags,nmin,nmax,&rest);
417: if (jac->reset) {
418: jac->head->is = rest;
419: jac->head->next->is = zerodiags;
420: } else {
421: PCFieldSplitSetIS(pc,"0",rest);
422: PCFieldSplitSetIS(pc,"1",zerodiags);
423: }
424: ISDestroy(&zerodiags);
425: ISDestroy(&rest);
426: } else if (coupling) {
427: IS coupling,rest;
428: PetscInt nmin,nmax;
430: MatGetOwnershipRange(pc->mat,&nmin,&nmax);
431: MatFindOffBlockDiagonalEntries(pc->mat,&coupling);
432: ISCreateStride(PetscObjectComm((PetscObject)pc->mat),nmax-nmin,nmin,1,&rest);
433: if (jac->reset) {
434: jac->head->is = coupling;
435: jac->head->next->is = rest;
436: } else {
437: PCFieldSplitSetIS(pc,"0",coupling);
438: PCFieldSplitSetIS(pc,"1",rest);
439: }
440: ISDestroy(&coupling);
441: ISDestroy(&rest);
442: } else {
443: if (jac->reset) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cases not yet handled when PCReset() was used");
444: PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&fieldsplit_default,NULL);
445: if (!fieldsplit_default) {
446: /* Allow user to set fields from command line, if bs was known at the time of PCSetFromOptions_FieldSplit()
447: then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */
448: PCFieldSplitSetRuntimeSplits_Private(pc);
449: if (jac->splitdefined) {PetscInfo(pc,"Splits defined using the options database\n");}
450: }
451: if (fieldsplit_default || !jac->splitdefined) {
452: PetscInfo(pc,"Using default splitting of fields\n");
453: for (i=0; i<jac->bs; i++) {
454: char splitname[8];
455: PetscSNPrintf(splitname,sizeof(splitname),"%D",i);
456: PCFieldSplitSetFields(pc,splitname,1,&i,&i);
457: }
458: jac->defaultsplit = PETSC_TRUE;
459: }
460: }
461: }
462: } else if (jac->nsplits == 1) {
463: if (ilink->is) {
464: IS is2;
465: PetscInt nmin,nmax;
467: MatGetOwnershipRange(pc->mat,&nmin,&nmax);
468: ISComplement(ilink->is,nmin,nmax,&is2);
469: PCFieldSplitSetIS(pc,"1",is2);
470: ISDestroy(&is2);
471: } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Must provide at least two sets of fields to PCFieldSplit()");
472: }
475: if (jac->nsplits < 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Unhandled case, must have at least two fields, not %d", jac->nsplits);
476: return(0);
477: }
479: PETSC_EXTERN PetscErrorCode PetscOptionsFindPairPrefix_Private(const char pre[], const char name[], char *value[], PetscBool *flg);
483: static PetscErrorCode PCSetUp_FieldSplit(PC pc)
484: {
485: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
486: PetscErrorCode ierr;
487: PC_FieldSplitLink ilink;
488: PetscInt i,nsplit;
489: PetscBool sorted, sorted_col;
492: PCFieldSplitSetDefaults(pc);
493: nsplit = jac->nsplits;
494: ilink = jac->head;
496: /* get the matrices for each split */
497: if (!jac->issetup) {
498: PetscInt rstart,rend,nslots,bs;
500: jac->issetup = PETSC_TRUE;
502: /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */
503: if (jac->defaultsplit || !ilink->is) {
504: if (jac->bs <= 0) jac->bs = nsplit;
505: }
506: bs = jac->bs;
507: MatGetOwnershipRange(pc->pmat,&rstart,&rend);
508: nslots = (rend - rstart)/bs;
509: for (i=0; i<nsplit; i++) {
510: if (jac->defaultsplit) {
511: ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+i,nsplit,&ilink->is);
512: ISDuplicate(ilink->is,&ilink->is_col);
513: } else if (!ilink->is) {
514: if (ilink->nfields > 1) {
515: PetscInt *ii,*jj,j,k,nfields = ilink->nfields,*fields = ilink->fields,*fields_col = ilink->fields_col;
516: PetscMalloc1(ilink->nfields*nslots,&ii);
517: PetscMalloc1(ilink->nfields*nslots,&jj);
518: for (j=0; j<nslots; j++) {
519: for (k=0; k<nfields; k++) {
520: ii[nfields*j + k] = rstart + bs*j + fields[k];
521: jj[nfields*j + k] = rstart + bs*j + fields_col[k];
522: }
523: }
524: ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,ii,PETSC_OWN_POINTER,&ilink->is);
525: ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,jj,PETSC_OWN_POINTER,&ilink->is_col);
526: } else {
527: ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields[0],bs,&ilink->is);
528: ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields_col[0],bs,&ilink->is_col);
529: }
530: }
531: ISSorted(ilink->is,&sorted);
532: if (ilink->is_col) { ISSorted(ilink->is_col,&sorted_col); }
533: if (!sorted || !sorted_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Fields must be sorted when creating split");
534: ilink = ilink->next;
535: }
536: }
538: ilink = jac->head;
539: if (!jac->pmat) {
540: Vec xtmp;
542: MatGetVecs(pc->pmat,&xtmp,NULL);
543: PetscMalloc1(nsplit,&jac->pmat);
544: PetscMalloc2(nsplit,&jac->x,nsplit,&jac->y);
545: for (i=0; i<nsplit; i++) {
546: MatNullSpace sp;
548: /* Check for preconditioning matrix attached to IS */
549: PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &jac->pmat[i]);
550: if (jac->pmat[i]) {
551: PetscObjectReference((PetscObject) jac->pmat[i]);
552: if (jac->type == PC_COMPOSITE_SCHUR) {
553: jac->schur_user = jac->pmat[i];
555: PetscObjectReference((PetscObject) jac->schur_user);
556: }
557: } else {
558: const char *prefix;
559: MatGetSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->pmat[i]);
560: KSPGetOptionsPrefix(ilink->ksp,&prefix);
561: MatSetOptionsPrefix(jac->pmat[i],prefix);
562: MatViewFromOptions(jac->pmat[i],NULL,"-mat_view");
563: }
564: /* create work vectors for each split */
565: MatGetVecs(jac->pmat[i],&jac->x[i],&jac->y[i]);
566: ilink->x = jac->x[i]; ilink->y = jac->y[i]; ilink->z = NULL;
567: /* compute scatter contexts needed by multiplicative versions and non-default splits */
568: VecScatterCreate(xtmp,ilink->is,jac->x[i],NULL,&ilink->sctx);
569: /* Check for null space attached to IS */
570: PetscObjectQuery((PetscObject) ilink->is, "nullspace", (PetscObject*) &sp);
571: if (sp) {
572: MatSetNullSpace(jac->pmat[i], sp);
573: }
574: PetscObjectQuery((PetscObject) ilink->is, "nearnullspace", (PetscObject*) &sp);
575: if (sp) {
576: MatSetNearNullSpace(jac->pmat[i], sp);
577: }
578: ilink = ilink->next;
579: }
580: VecDestroy(&xtmp);
581: } else {
582: for (i=0; i<nsplit; i++) {
583: Mat pmat;
585: /* Check for preconditioning matrix attached to IS */
586: PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &pmat);
587: if (!pmat) {
588: MatGetSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_REUSE_MATRIX,&jac->pmat[i]);
589: }
590: ilink = ilink->next;
591: }
592: }
593: if (jac->diag_use_amat) {
594: ilink = jac->head;
595: if (!jac->mat) {
596: PetscMalloc1(nsplit,&jac->mat);
597: for (i=0; i<nsplit; i++) {
598: MatGetSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->mat[i]);
599: ilink = ilink->next;
600: }
601: } else {
602: for (i=0; i<nsplit; i++) {
603: if (jac->mat[i]) {MatGetSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_REUSE_MATRIX,&jac->mat[i]);}
604: ilink = ilink->next;
605: }
606: }
607: } else {
608: jac->mat = jac->pmat;
609: }
611: if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR) {
612: /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */
613: /* FIXME: Can/should we reuse jac->mat whenever (jac->diag_use_amat) is true? */
614: ilink = jac->head;
615: if (!jac->Afield) {
616: PetscMalloc1(nsplit,&jac->Afield);
617: for (i=0; i<nsplit; i++) {
618: MatGetSubMatrix(pc->mat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);
619: ilink = ilink->next;
620: }
621: } else {
622: for (i=0; i<nsplit; i++) {
623: MatGetSubMatrix(pc->mat,ilink->is,NULL,MAT_REUSE_MATRIX,&jac->Afield[i]);
624: ilink = ilink->next;
625: }
626: }
627: }
629: if (jac->type == PC_COMPOSITE_SCHUR) {
630: IS ccis;
631: PetscInt rstart,rend;
632: char lscname[256];
633: PetscObject LSC_L;
635: if (nsplit != 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"To use Schur complement preconditioner you must have exactly 2 fields");
637: /* When extracting off-diagonal submatrices, we take complements from this range */
638: MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);
640: /* need to handle case when one is resetting up the preconditioner */
641: if (jac->schur) {
642: KSP kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper;
644: MatSchurComplementGetKSP(jac->schur, &kspInner);
645: ilink = jac->head;
646: ISComplement(ilink->is_col,rstart,rend,&ccis);
647: if (jac->offdiag_use_amat) {
648: MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);
649: } else {
650: MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);
651: }
652: ISDestroy(&ccis);
653: ilink = ilink->next;
654: ISComplement(ilink->is_col,rstart,rend,&ccis);
655: if (jac->offdiag_use_amat) {
656: MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);
657: } else {
658: MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);
659: }
660: ISDestroy(&ccis);
661: MatSchurComplementUpdateSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);
662: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
663: MatDestroy(&jac->schurp);
664: MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);
665: }
666: if (kspA != kspInner) {
667: KSPSetOperators(kspA,jac->mat[0],jac->pmat[0]);
668: }
669: if (kspUpper != kspA) {
670: KSPSetOperators(kspUpper,jac->mat[0],jac->pmat[0]);
671: }
672: KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));
673: } else {
674: const char *Dprefix;
675: char schurprefix[256], schurmatprefix[256];
676: char schurtestoption[256];
677: MatNullSpace sp;
678: PetscBool flg;
680: /* extract the A01 and A10 matrices */
681: ilink = jac->head;
682: ISComplement(ilink->is_col,rstart,rend,&ccis);
683: if (jac->offdiag_use_amat) {
684: MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
685: } else {
686: MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
687: }
688: ISDestroy(&ccis);
689: ilink = ilink->next;
690: ISComplement(ilink->is_col,rstart,rend,&ccis);
691: if (jac->offdiag_use_amat) {
692: MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
693: } else {
694: MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
695: }
696: ISDestroy(&ccis);
698: /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */
699: MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur);
700: MatSetType(jac->schur,MATSCHURCOMPLEMENT);
701: MatSchurComplementSetSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);
702: PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
703: /* 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? */
704: MatSetOptionsPrefix(jac->schur,schurmatprefix);
705: MatSetFromOptions(jac->schur);
706: MatGetNullSpace(jac->pmat[1], &sp);
707: if (sp) {
708: MatSetNullSpace(jac->schur, sp);
709: }
711: PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);
712: PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
713: if (flg) {
714: DM dmInner;
715: KSP kspInner;
717: MatSchurComplementGetKSP(jac->schur, &kspInner);
718: PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
719: /* Indent this deeper to emphasize the "inner" nature of this solver. */
720: PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject) pc, 2);
721: KSPSetOptionsPrefix(kspInner, schurprefix);
723: /* Set DM for new solver */
724: KSPGetDM(jac->head->ksp, &dmInner);
725: KSPSetDM(kspInner, dmInner);
726: KSPSetDMActive(kspInner, PETSC_FALSE);
727: } else {
728: /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or
729: * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact,
730: * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for
731: * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make
732: * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used
733: * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */
734: KSPSetType(jac->head->ksp,KSPGMRES);
735: MatSchurComplementSetKSP(jac->schur,jac->head->ksp);
736: }
737: KSPSetOperators(jac->head->ksp,jac->mat[0],jac->pmat[0]);
738: KSPSetFromOptions(jac->head->ksp);
739: MatSetFromOptions(jac->schur);
741: PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);
742: PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
743: if (flg) {
744: DM dmInner;
746: PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
747: KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper);
748: KSPSetOptionsPrefix(jac->kspupper, schurprefix);
749: KSPGetDM(jac->head->ksp, &dmInner);
750: KSPSetDM(jac->kspupper, dmInner);
751: KSPSetDMActive(jac->kspupper, PETSC_FALSE);
752: KSPSetFromOptions(jac->kspupper);
753: KSPSetOperators(jac->kspupper,jac->mat[0],jac->pmat[0]);
754: VecDuplicate(jac->head->x, &jac->head->z);
755: } else {
756: jac->kspupper = jac->head->ksp;
757: PetscObjectReference((PetscObject) jac->head->ksp);
758: }
760: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
761: MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);
762: }
763: KSPCreate(PetscObjectComm((PetscObject)pc),&jac->kspschur);
764: PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);
765: PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);
766: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) {
767: PC pcschur;
768: KSPGetPC(jac->kspschur,&pcschur);
769: PCSetType(pcschur,PCNONE);
770: /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */
771: } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) {
772: MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user);
773: }
774: KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));
775: KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix);
776: KSPSetOptionsPrefix(jac->kspschur, Dprefix);
777: /* propogate DM */
778: {
779: DM sdm;
780: KSPGetDM(jac->head->next->ksp, &sdm);
781: if (sdm) {
782: KSPSetDM(jac->kspschur, sdm);
783: KSPSetDMActive(jac->kspschur, PETSC_FALSE);
784: }
785: }
786: /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
787: /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */
788: KSPSetFromOptions(jac->kspschur);
789: }
791: /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */
792: PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_L",ilink->splitname);
793: PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);
794: if (!LSC_L) {PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);}
795: if (LSC_L) {PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L);}
796: PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_Lp",ilink->splitname);
797: PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);
798: if (!LSC_L) {PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);}
799: if (LSC_L) {PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L);}
800: } else {
801: /* set up the individual splits' PCs */
802: i = 0;
803: ilink = jac->head;
804: while (ilink) {
805: KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i]);
806: /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
807: if (!jac->suboptionsset) {KSPSetFromOptions(ilink->ksp);}
808: i++;
809: ilink = ilink->next;
810: }
811: }
813: jac->suboptionsset = PETSC_TRUE;
814: return(0);
815: }
817: #define FieldSplitSplitSolveAdd(ilink,xx,yy) \
818: (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
819: VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
820: KSPSolve(ilink->ksp,ilink->x,ilink->y) || \
821: VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \
822: VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE))
826: static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y)
827: {
828: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
829: PetscErrorCode ierr;
830: PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
831: KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;
834: switch (jac->schurfactorization) {
835: case PC_FIELDSPLIT_SCHUR_FACT_DIAG:
836: /* [A00 0; 0 -S], positive definite, suitable for MINRES */
837: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
838: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
839: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
840: KSPSolve(kspA,ilinkA->x,ilinkA->y);
841: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
842: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
843: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
844: VecScale(ilinkD->y,-1.);
845: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
846: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
847: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
848: break;
849: case PC_FIELDSPLIT_SCHUR_FACT_LOWER:
850: /* [A00 0; A10 S], suitable for left preconditioning */
851: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
852: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
853: KSPSolve(kspA,ilinkA->x,ilinkA->y);
854: MatMult(jac->C,ilinkA->y,ilinkD->x);
855: VecScale(ilinkD->x,-1.);
856: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
857: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
858: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
859: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
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_UPPER:
865: /* [A00 A01; 0 S], suitable for right preconditioning */
866: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
867: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
868: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
869: MatMult(jac->B,ilinkD->y,ilinkA->x);
870: VecScale(ilinkA->x,-1.);
871: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);
872: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
873: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);
874: KSPSolve(kspA,ilinkA->x,ilinkA->y);
875: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
876: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
877: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
878: break;
879: case PC_FIELDSPLIT_SCHUR_FACT_FULL:
880: /* [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 */
881: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
882: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
883: KSPSolve(kspLower,ilinkA->x,ilinkA->y);
884: MatMult(jac->C,ilinkA->y,ilinkD->x);
885: VecScale(ilinkD->x,-1.0);
886: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
887: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
889: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
890: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
891: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
893: if (kspUpper == kspA) {
894: MatMult(jac->B,ilinkD->y,ilinkA->y);
895: VecAXPY(ilinkA->x,-1.0,ilinkA->y);
896: KSPSolve(kspA,ilinkA->x,ilinkA->y);
897: } else {
898: KSPSolve(kspA,ilinkA->x,ilinkA->y);
899: MatMult(jac->B,ilinkD->y,ilinkA->x);
900: KSPSolve(kspUpper,ilinkA->x,ilinkA->z);
901: VecAXPY(ilinkA->y,-1.0,ilinkA->z);
902: }
903: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
904: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
905: }
906: return(0);
907: }
911: static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y)
912: {
913: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
914: PetscErrorCode ierr;
915: PC_FieldSplitLink ilink = jac->head;
916: PetscInt cnt,bs;
919: if (jac->type == PC_COMPOSITE_ADDITIVE) {
920: if (jac->defaultsplit) {
921: VecGetBlockSize(x,&bs);
922: 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);
923: VecGetBlockSize(y,&bs);
924: 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);
925: VecStrideGatherAll(x,jac->x,INSERT_VALUES);
926: while (ilink) {
927: KSPSolve(ilink->ksp,ilink->x,ilink->y);
928: ilink = ilink->next;
929: }
930: VecStrideScatterAll(jac->y,y,INSERT_VALUES);
931: } else {
932: VecSet(y,0.0);
933: while (ilink) {
934: FieldSplitSplitSolveAdd(ilink,x,y);
935: ilink = ilink->next;
936: }
937: }
938: } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
939: if (!jac->w1) {
940: VecDuplicate(x,&jac->w1);
941: VecDuplicate(x,&jac->w2);
942: }
943: VecSet(y,0.0);
944: FieldSplitSplitSolveAdd(ilink,x,y);
945: cnt = 1;
946: while (ilink->next) {
947: ilink = ilink->next;
948: /* compute the residual only over the part of the vector needed */
949: MatMult(jac->Afield[cnt++],y,ilink->x);
950: VecScale(ilink->x,-1.0);
951: VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
952: VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
953: KSPSolve(ilink->ksp,ilink->x,ilink->y);
954: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
955: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
956: }
957: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
958: cnt -= 2;
959: while (ilink->previous) {
960: ilink = ilink->previous;
961: /* compute the residual only over the part of the vector needed */
962: MatMult(jac->Afield[cnt--],y,ilink->x);
963: VecScale(ilink->x,-1.0);
964: VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
965: VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
966: KSPSolve(ilink->ksp,ilink->x,ilink->y);
967: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
968: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
969: }
970: }
971: } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type);
972: return(0);
973: }
975: #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \
976: (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
977: VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
978: KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \
979: VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \
980: VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE))
984: static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y)
985: {
986: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
987: PetscErrorCode ierr;
988: PC_FieldSplitLink ilink = jac->head;
989: PetscInt bs;
992: if (jac->type == PC_COMPOSITE_ADDITIVE) {
993: if (jac->defaultsplit) {
994: VecGetBlockSize(x,&bs);
995: 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);
996: VecGetBlockSize(y,&bs);
997: 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);
998: VecStrideGatherAll(x,jac->x,INSERT_VALUES);
999: while (ilink) {
1000: KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);
1001: ilink = ilink->next;
1002: }
1003: VecStrideScatterAll(jac->y,y,INSERT_VALUES);
1004: } else {
1005: VecSet(y,0.0);
1006: while (ilink) {
1007: FieldSplitSplitSolveAddTranspose(ilink,x,y);
1008: ilink = ilink->next;
1009: }
1010: }
1011: } else {
1012: if (!jac->w1) {
1013: VecDuplicate(x,&jac->w1);
1014: VecDuplicate(x,&jac->w2);
1015: }
1016: VecSet(y,0.0);
1017: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1018: FieldSplitSplitSolveAddTranspose(ilink,x,y);
1019: while (ilink->next) {
1020: ilink = ilink->next;
1021: MatMultTranspose(pc->mat,y,jac->w1);
1022: VecWAXPY(jac->w2,-1.0,jac->w1,x);
1023: FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1024: }
1025: while (ilink->previous) {
1026: ilink = ilink->previous;
1027: MatMultTranspose(pc->mat,y,jac->w1);
1028: VecWAXPY(jac->w2,-1.0,jac->w1,x);
1029: FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1030: }
1031: } else {
1032: while (ilink->next) { /* get to last entry in linked list */
1033: ilink = ilink->next;
1034: }
1035: FieldSplitSplitSolveAddTranspose(ilink,x,y);
1036: while (ilink->previous) {
1037: ilink = ilink->previous;
1038: MatMultTranspose(pc->mat,y,jac->w1);
1039: VecWAXPY(jac->w2,-1.0,jac->w1,x);
1040: FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1041: }
1042: }
1043: }
1044: return(0);
1045: }
1049: static PetscErrorCode PCReset_FieldSplit(PC pc)
1050: {
1051: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1052: PetscErrorCode ierr;
1053: PC_FieldSplitLink ilink = jac->head,next;
1056: while (ilink) {
1057: KSPReset(ilink->ksp);
1058: VecDestroy(&ilink->x);
1059: VecDestroy(&ilink->y);
1060: VecDestroy(&ilink->z);
1061: VecScatterDestroy(&ilink->sctx);
1062: ISDestroy(&ilink->is);
1063: ISDestroy(&ilink->is_col);
1064: next = ilink->next;
1065: ilink = next;
1066: }
1067: PetscFree2(jac->x,jac->y);
1068: if (jac->mat && jac->mat != jac->pmat) {
1069: MatDestroyMatrices(jac->nsplits,&jac->mat);
1070: } else if (jac->mat) {
1071: jac->mat = NULL;
1072: }
1073: if (jac->pmat) {MatDestroyMatrices(jac->nsplits,&jac->pmat);}
1074: if (jac->Afield) {MatDestroyMatrices(jac->nsplits,&jac->Afield);}
1075: VecDestroy(&jac->w1);
1076: VecDestroy(&jac->w2);
1077: MatDestroy(&jac->schur);
1078: MatDestroy(&jac->schurp);
1079: MatDestroy(&jac->schur_user);
1080: KSPDestroy(&jac->kspschur);
1081: KSPDestroy(&jac->kspupper);
1082: MatDestroy(&jac->B);
1083: MatDestroy(&jac->C);
1084: jac->reset = PETSC_TRUE;
1085: return(0);
1086: }
1090: static PetscErrorCode PCDestroy_FieldSplit(PC pc)
1091: {
1092: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1093: PetscErrorCode ierr;
1094: PC_FieldSplitLink ilink = jac->head,next;
1097: PCReset_FieldSplit(pc);
1098: while (ilink) {
1099: KSPDestroy(&ilink->ksp);
1100: next = ilink->next;
1101: PetscFree(ilink->splitname);
1102: PetscFree(ilink->fields);
1103: PetscFree(ilink->fields_col);
1104: PetscFree(ilink);
1105: ilink = next;
1106: }
1107: PetscFree2(jac->x,jac->y);
1108: PetscFree(pc->data);
1109: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL);
1110: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL);
1111: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL);
1112: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL);
1113: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL);
1114: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL);
1115: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL);
1116: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL);
1117: return(0);
1118: }
1122: static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc)
1123: {
1124: PetscErrorCode ierr;
1125: PetscInt bs;
1126: PetscBool flg,stokes = PETSC_FALSE;
1127: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1128: PCCompositeType ctype;
1131: PetscOptionsHead("FieldSplit options");
1132: PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL);
1133: PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);
1134: if (flg) {
1135: PCFieldSplitSetBlockSize(pc,bs);
1136: }
1138: PetscOptionsBool("-pc_fieldsplit_diag_use_amat","Use Amat (not Pmat) to extract diagonal fieldsplit blocks", "PCFieldSplitSetDiagUseAmat",pc->useAmat,&jac->diag_use_amat,NULL);
1139: PetscOptionsBool("-pc_fieldsplit_off_diag_use_amat","Use Amat (not Pmat) to extract off-diagonal fieldsplit blocks", "PCFieldSplitSetOffDiagUseAmat",pc->useAmat,&jac->offdiag_use_amat,NULL);
1140: /* FIXME: No programmatic equivalent to the following. */
1141: PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);
1142: if (stokes) {
1143: PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);
1144: jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF;
1145: }
1147: PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);
1148: if (flg) {
1149: PCFieldSplitSetType(pc,ctype);
1150: }
1151: /* Only setup fields once */
1152: if ((jac->bs > 0) && (jac->nsplits == 0)) {
1153: /* only allow user to set fields from command line if bs is already known.
1154: otherwise user can set them in PCFieldSplitSetDefaults() */
1155: PCFieldSplitSetRuntimeSplits_Private(pc);
1156: if (jac->splitdefined) {PetscInfo(pc,"Splits defined using the options database\n");}
1157: }
1158: if (jac->type == PC_COMPOSITE_SCHUR) {
1159: PetscOptionsGetEnum(((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);
1160: if (flg) {PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");}
1161: 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);
1162: PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSetSchurPre",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL);
1163: }
1164: PetscOptionsTail();
1165: return(0);
1166: }
1168: /*------------------------------------------------------------------------------------*/
1172: static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1173: {
1174: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1175: PetscErrorCode ierr;
1176: PC_FieldSplitLink ilink,next = jac->head;
1177: char prefix[128];
1178: PetscInt i;
1181: if (jac->splitdefined) {
1182: PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1183: return(0);
1184: }
1185: for (i=0; i<n; i++) {
1186: if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs);
1187: if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]);
1188: }
1189: PetscNew(&ilink);
1190: if (splitname) {
1191: PetscStrallocpy(splitname,&ilink->splitname);
1192: } else {
1193: PetscMalloc1(3,&ilink->splitname);
1194: PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);
1195: }
1196: PetscMalloc1(n,&ilink->fields);
1197: PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));
1198: PetscMalloc1(n,&ilink->fields_col);
1199: PetscMemcpy(ilink->fields_col,fields_col,n*sizeof(PetscInt));
1201: ilink->nfields = n;
1202: ilink->next = NULL;
1203: KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1204: PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1205: KSPSetType(ilink->ksp,KSPPREONLY);
1206: PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);
1208: PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);
1209: KSPSetOptionsPrefix(ilink->ksp,prefix);
1211: if (!next) {
1212: jac->head = ilink;
1213: ilink->previous = NULL;
1214: } else {
1215: while (next->next) {
1216: next = next->next;
1217: }
1218: next->next = ilink;
1219: ilink->previous = next;
1220: }
1221: jac->nsplits++;
1222: return(0);
1223: }
1227: static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp)
1228: {
1229: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1233: PetscMalloc1(jac->nsplits,subksp);
1234: MatSchurComplementGetKSP(jac->schur,*subksp);
1236: (*subksp)[1] = jac->kspschur;
1237: if (n) *n = jac->nsplits;
1238: return(0);
1239: }
1243: static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp)
1244: {
1245: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1246: PetscErrorCode ierr;
1247: PetscInt cnt = 0;
1248: PC_FieldSplitLink ilink = jac->head;
1251: PetscMalloc1(jac->nsplits,subksp);
1252: while (ilink) {
1253: (*subksp)[cnt++] = ilink->ksp;
1254: ilink = ilink->next;
1255: }
1256: 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);
1257: if (n) *n = jac->nsplits;
1258: return(0);
1259: }
1263: static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is)
1264: {
1265: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1266: PetscErrorCode ierr;
1267: PC_FieldSplitLink ilink, next = jac->head;
1268: char prefix[128];
1271: if (jac->splitdefined) {
1272: PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1273: return(0);
1274: }
1275: PetscNew(&ilink);
1276: if (splitname) {
1277: PetscStrallocpy(splitname,&ilink->splitname);
1278: } else {
1279: PetscMalloc1(8,&ilink->splitname);
1280: PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);
1281: }
1282: PetscObjectReference((PetscObject)is);
1283: ISDestroy(&ilink->is);
1284: ilink->is = is;
1285: PetscObjectReference((PetscObject)is);
1286: ISDestroy(&ilink->is_col);
1287: ilink->is_col = is;
1288: ilink->next = NULL;
1289: KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1290: PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1291: KSPSetType(ilink->ksp,KSPPREONLY);
1292: PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);
1294: PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);
1295: KSPSetOptionsPrefix(ilink->ksp,prefix);
1297: if (!next) {
1298: jac->head = ilink;
1299: ilink->previous = NULL;
1300: } else {
1301: while (next->next) {
1302: next = next->next;
1303: }
1304: next->next = ilink;
1305: ilink->previous = next;
1306: }
1307: jac->nsplits++;
1308: return(0);
1309: }
1313: /*@
1314: PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner
1316: Logically Collective on PC
1318: Input Parameters:
1319: + pc - the preconditioner context
1320: . splitname - name of this split, if NULL the number of the split is used
1321: . n - the number of fields in this split
1322: - fields - the fields in this split
1324: Level: intermediate
1326: Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field.
1328: The PCFieldSplitSetFields() is for defining fields as strided blocks. For example, if the block
1329: 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
1330: 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x....
1331: where the numbered entries indicate what is in the field.
1333: This function is called once per split (it creates a new split each time). Solve options
1334: for this split will be available under the prefix -fieldsplit_SPLITNAME_.
1336: Developer Note: This routine does not actually create the IS representing the split, that is delayed
1337: until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be
1338: available when this routine is called.
1340: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS()
1342: @*/
1343: PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1344: {
1350: if (n < 1) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname);
1352: PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt*,const PetscInt*),(pc,splitname,n,fields,fields_col));
1353: return(0);
1354: }
1358: /*@
1359: PCFieldSplitSetDiagUseAmat - set flag indicating whether to extract diagonal blocks from Amat (rather than Pmat)
1361: Logically Collective on PC
1363: Input Parameters:
1364: + pc - the preconditioner object
1365: - flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from
1367: Options Database:
1368: . -pc_fieldsplit_diag_use_amat
1370: Level: intermediate
1372: .seealso: PCFieldSplitGetDiagUseAmat(), PCFieldSplitSetOffDiagUseAmat(), PCFIELDSPLIT
1374: @*/
1375: PetscErrorCode PCFieldSplitSetDiagUseAmat(PC pc,PetscBool flg)
1376: {
1377: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1378: PetscBool isfs;
1383: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1384: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1385: jac->diag_use_amat = flg;
1386: return(0);
1387: }
1391: /*@
1392: PCFieldSplitGetDiagUseAmat - get the flag indicating whether to extract diagonal blocks from Amat (rather than Pmat)
1394: Logically Collective on PC
1396: Input Parameters:
1397: . pc - the preconditioner object
1399: Output Parameters:
1400: . flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from
1403: Level: intermediate
1405: .seealso: PCFieldSplitSetDiagUseAmat(), PCFieldSplitGetOffDiagUseAmat(), PCFIELDSPLIT
1407: @*/
1408: PetscErrorCode PCFieldSplitGetDiagUseAmat(PC pc,PetscBool *flg)
1409: {
1410: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1411: PetscBool isfs;
1417: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1418: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1419: *flg = jac->diag_use_amat;
1420: return(0);
1421: }
1425: /*@
1426: PCFieldSplitSetOffDiagUseAmat - set flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat)
1428: Logically Collective on PC
1430: Input Parameters:
1431: + pc - the preconditioner object
1432: - flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from
1434: Options Database:
1435: . -pc_fieldsplit_off_diag_use_amat
1437: Level: intermediate
1439: .seealso: PCFieldSplitGetOffDiagUseAmat(), PCFieldSplitSetDiagUseAmat(), PCFIELDSPLIT
1441: @*/
1442: PetscErrorCode PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg)
1443: {
1444: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1445: PetscBool isfs;
1450: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1451: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1452: jac->offdiag_use_amat = flg;
1453: return(0);
1454: }
1458: /*@
1459: PCFieldSplitGetOffDiagUseAmat - get the flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat)
1461: Logically Collective on PC
1463: Input Parameters:
1464: . pc - the preconditioner object
1466: Output Parameters:
1467: . flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from
1470: Level: intermediate
1472: .seealso: PCFieldSplitSetOffDiagUseAmat(), PCFieldSplitGetDiagUseAmat(), PCFIELDSPLIT
1474: @*/
1475: PetscErrorCode PCFieldSplitGetOffDiagUseAmat(PC pc,PetscBool *flg)
1476: {
1477: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1478: PetscBool isfs;
1484: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1485: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1486: *flg = jac->offdiag_use_amat;
1487: return(0);
1488: }
1494: /*@
1495: PCFieldSplitSetIS - Sets the exact elements for field
1497: Logically Collective on PC
1499: Input Parameters:
1500: + pc - the preconditioner context
1501: . splitname - name of this split, if NULL the number of the split is used
1502: - is - the index set that defines the vector elements in this field
1505: Notes:
1506: Use PCFieldSplitSetFields(), for fields defined by strided types.
1508: This function is called once per split (it creates a new split each time). Solve options
1509: for this split will be available under the prefix -fieldsplit_SPLITNAME_.
1511: Level: intermediate
1513: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize()
1515: @*/
1516: PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is)
1517: {
1524: PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));
1525: return(0);
1526: }
1530: /*@
1531: PCFieldSplitGetIS - Retrieves the elements for a field as an IS
1533: Logically Collective on PC
1535: Input Parameters:
1536: + pc - the preconditioner context
1537: - splitname - name of this split
1539: Output Parameter:
1540: - is - the index set that defines the vector elements in this field, or NULL if the field is not found
1542: Level: intermediate
1544: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS()
1546: @*/
1547: PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is)
1548: {
1555: {
1556: PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;
1557: PC_FieldSplitLink ilink = jac->head;
1558: PetscBool found;
1560: *is = NULL;
1561: while (ilink) {
1562: PetscStrcmp(ilink->splitname, splitname, &found);
1563: if (found) {
1564: *is = ilink->is;
1565: break;
1566: }
1567: ilink = ilink->next;
1568: }
1569: }
1570: return(0);
1571: }
1575: /*@
1576: PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the
1577: fieldsplit preconditioner. If not set the matrix block size is used.
1579: Logically Collective on PC
1581: Input Parameters:
1582: + pc - the preconditioner context
1583: - bs - the block size
1585: Level: intermediate
1587: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields()
1589: @*/
1590: PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs)
1591: {
1597: PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));
1598: return(0);
1599: }
1603: /*@C
1604: PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits
1606: Collective on KSP
1608: Input Parameter:
1609: . pc - the preconditioner context
1611: Output Parameters:
1612: + n - the number of splits
1613: - pc - the array of KSP contexts
1615: Note:
1616: After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user
1617: (not the KSP just the array that contains them).
1619: You must call KSPSetUp() before calling PCFieldSplitGetSubKSP().
1621: Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs.
1622: You can call PCFieldSplitGetSubKSP(pc,n,NULL_OBJECT,ierr) to determine how large the
1623: KSP array must be.
1626: Level: advanced
1628: .seealso: PCFIELDSPLIT
1629: @*/
1630: PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
1631: {
1637: PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));
1638: return(0);
1639: }
1643: /*@
1644: PCFieldSplitSetSchurPre - Indicates if the Schur complement is preconditioned by a preconditioner constructed by the
1645: A11 matrix. Otherwise no preconditioner is used.
1647: Collective on PC
1649: Input Parameters:
1650: + pc - the preconditioner context
1651: . 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
1652: - userpre - matrix to use for preconditioning, or NULL
1654: Options Database:
1655: . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> default is a11
1657: Notes:
1658: $ If ptype is
1659: $ user then the preconditioner for the Schur complement is generated by the provided matrix (pre argument
1660: $ to this function).
1661: $ a11 then the preconditioner for the Schur complement is generated by the block diagonal part of the preconditioner
1662: $ matrix associated with the Schur complement (i.e. A11), not he Schur complement matrix
1663: $ full then the preconditioner uses the exact Schur complement (this is expensive)
1664: $ self the preconditioner for the Schur complement is generated from the Schur complement matrix itself:
1665: $ The only preconditioner that currently works directly with the Schur complement matrix object is the PCLSC
1666: $ preconditioner
1667: $ selfp then the preconditioning matrix is an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01
1668: $ This is only a good preconditioner when diag(A00) is a good preconditioner for A00.
1670: When solving a saddle point problem, where the A11 block is identically zero, using a11 as the ptype only makes sense
1671: with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and
1672: -fieldsplit_1_pc_type lsc which uses the least squares commutator to compute a preconditioner for the Schur complement.
1674: Level: intermediate
1676: .seealso: PCFieldSplitGetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC
1678: @*/
1679: PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
1680: {
1685: PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));
1686: return(0);
1687: }
1688: PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) {return PCFieldSplitSetSchurPre(pc,ptype,pre);} /* Deprecated name */
1692: /*@
1693: PCFieldSplitGetSchurPre - For Schur complement fieldsplit, determine how the Schur complement will be
1694: preconditioned. See PCFieldSplitSetSchurPre() for details.
1696: Logically Collective on PC
1698: Input Parameters:
1699: . pc - the preconditioner context
1701: Output Parameters:
1702: + ptype - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11, PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_PRE_USER
1703: - userpre - matrix to use for preconditioning (with PC_FIELDSPLIT_PRE_USER), or NULL
1705: Level: intermediate
1707: .seealso: PCFieldSplitSetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC
1709: @*/
1710: PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
1711: {
1716: PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre));
1717: return(0);
1718: }
1722: /*@
1723: PCFieldSplitSchurGetS - extract the MatSchurComplement object used by this PC in case it needs to be configured separately
1725: Not collective
1727: Input Parameter:
1728: . pc - the preconditioner context
1730: Output Parameter:
1731: . S - the Schur complement matrix
1733: Notes:
1734: This matrix should not be destroyed using MatDestroy(); rather, use PCFieldSplitSchurRestoreS().
1736: Level: advanced
1738: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurRestoreS()
1740: @*/
1741: PetscErrorCode PCFieldSplitSchurGetS(PC pc,Mat *S)
1742: {
1744: const char* t;
1745: PetscBool isfs;
1746: PC_FieldSplit *jac;
1750: PetscObjectGetType((PetscObject)pc,&t);
1751: PetscStrcmp(t,PCFIELDSPLIT,&isfs);
1752: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
1753: jac = (PC_FieldSplit*)pc->data;
1754: if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
1755: if (S) *S = jac->schur;
1756: return(0);
1757: }
1761: /*@
1762: PCFieldSplitSchurRestoreS - restores the MatSchurComplement object used by this PC
1764: Not collective
1766: Input Parameters:
1767: + pc - the preconditioner context
1768: . S - the Schur complement matrix
1770: Level: advanced
1772: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurGetS()
1774: @*/
1775: PetscErrorCode PCFieldSplitSchurRestoreS(PC pc,Mat *S)
1776: {
1778: const char* t;
1779: PetscBool isfs;
1780: PC_FieldSplit *jac;
1784: PetscObjectGetType((PetscObject)pc,&t);
1785: PetscStrcmp(t,PCFIELDSPLIT,&isfs);
1786: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
1787: jac = (PC_FieldSplit*)pc->data;
1788: if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
1789: if (!S || *S != jac->schur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten");
1790: return(0);
1791: }
1796: static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
1797: {
1798: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1802: jac->schurpre = ptype;
1803: if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) {
1804: MatDestroy(&jac->schur_user);
1805: jac->schur_user = pre;
1806: PetscObjectReference((PetscObject)jac->schur_user);
1807: }
1808: return(0);
1809: }
1813: static PetscErrorCode PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
1814: {
1815: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1818: *ptype = jac->schurpre;
1819: *pre = jac->schur_user;
1820: return(0);
1821: }
1825: /*@
1826: PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain
1828: Collective on PC
1830: Input Parameters:
1831: + pc - the preconditioner context
1832: - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default
1834: Options Database:
1835: . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full
1838: Level: intermediate
1840: Notes:
1841: The FULL factorization is
1843: $ (A B) = (1 0) (A 0) (1 Ainv*B)
1844: $ (C D) (C*Ainv 1) (0 S) (0 1 )
1846: 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,
1847: 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).
1849: 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
1850: 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
1851: 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,
1852: the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so KSPGMRES converges in at most 4 iterations.
1854: 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
1855: 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).
1857: References:
1858: Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000) pp. 1969-1972.
1860: Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001), pp. 1050-1051.
1862: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType
1863: @*/
1864: PetscErrorCode PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype)
1865: {
1870: PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));
1871: return(0);
1872: }
1876: static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype)
1877: {
1878: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1881: jac->schurfactorization = ftype;
1882: return(0);
1883: }
1887: /*@C
1888: PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement
1890: Collective on KSP
1892: Input Parameter:
1893: . pc - the preconditioner context
1895: Output Parameters:
1896: + A00 - the (0,0) block
1897: . A01 - the (0,1) block
1898: . A10 - the (1,0) block
1899: - A11 - the (1,1) block
1901: Level: advanced
1903: .seealso: PCFIELDSPLIT
1904: @*/
1905: PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11)
1906: {
1907: PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;
1911: if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach.");
1912: if (A00) *A00 = jac->pmat[0];
1913: if (A01) *A01 = jac->B;
1914: if (A10) *A10 = jac->C;
1915: if (A11) *A11 = jac->pmat[1];
1916: return(0);
1917: }
1921: static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type)
1922: {
1923: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1927: jac->type = type;
1928: if (type == PC_COMPOSITE_SCHUR) {
1929: pc->ops->apply = PCApply_FieldSplit_Schur;
1930: pc->ops->view = PCView_FieldSplit_Schur;
1932: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);
1933: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit);
1934: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit);
1935: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);
1937: } else {
1938: pc->ops->apply = PCApply_FieldSplit;
1939: pc->ops->view = PCView_FieldSplit;
1941: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
1942: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",0);
1943: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",0);
1944: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);
1945: }
1946: return(0);
1947: }
1951: static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs)
1952: {
1953: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1956: if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs);
1957: 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);
1958: jac->bs = bs;
1959: return(0);
1960: }
1964: /*@
1965: PCFieldSplitSetType - Sets the type of fieldsplit preconditioner.
1967: Collective on PC
1969: Input Parameter:
1970: . pc - the preconditioner context
1971: . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR
1973: Options Database Key:
1974: . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type
1976: Level: Intermediate
1978: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
1980: .seealso: PCCompositeSetType()
1982: @*/
1983: PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type)
1984: {
1989: PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));
1990: return(0);
1991: }
1995: /*@
1996: PCFieldSplitGetType - Gets the type of fieldsplit preconditioner.
1998: Not collective
2000: Input Parameter:
2001: . pc - the preconditioner context
2003: Output Parameter:
2004: . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR
2006: Level: Intermediate
2008: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
2009: .seealso: PCCompositeSetType()
2010: @*/
2011: PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type)
2012: {
2013: PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;
2018: *type = jac->type;
2019: return(0);
2020: }
2024: /*@
2025: PCFieldSplitSetDMSplits - Flags whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible.
2027: Logically Collective
2029: Input Parameters:
2030: + pc - the preconditioner context
2031: - flg - boolean indicating whether to use field splits defined by the DM
2033: Options Database Key:
2034: . -pc_fieldsplit_dm_splits
2036: Level: Intermediate
2038: .keywords: PC, DM, composite preconditioner, additive, multiplicative
2040: .seealso: PCFieldSplitGetDMSplits()
2042: @*/
2043: PetscErrorCode PCFieldSplitSetDMSplits(PC pc,PetscBool flg)
2044: {
2045: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2046: PetscBool isfs;
2052: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2053: if (isfs) {
2054: jac->dm_splits = flg;
2055: }
2056: return(0);
2057: }
2062: /*@
2063: PCFieldSplitGetDMSplits - Returns flag indicating whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible.
2065: Logically Collective
2067: Input Parameter:
2068: . pc - the preconditioner context
2070: Output Parameter:
2071: . flg - boolean indicating whether to use field splits defined by the DM
2073: Level: Intermediate
2075: .keywords: PC, DM, composite preconditioner, additive, multiplicative
2077: .seealso: PCFieldSplitSetDMSplits()
2079: @*/
2080: PetscErrorCode PCFieldSplitGetDMSplits(PC pc,PetscBool* flg)
2081: {
2082: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2083: PetscBool isfs;
2089: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2090: if (isfs) {
2091: if(flg) *flg = jac->dm_splits;
2092: }
2093: return(0);
2094: }
2096: /* -------------------------------------------------------------------------------------*/
2097: /*MC
2098: PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual
2099: fields or groups of fields. See the users manual section "Solving Block Matrices" for more details.
2101: To set options on the solvers for each block append -fieldsplit_ to all the PC
2102: options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1
2104: To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP()
2105: and set the options directly on the resulting KSP object
2107: Level: intermediate
2109: Options Database Keys:
2110: + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split
2111: . -pc_fieldsplit_default - automatically add any fields to additional splits that have not
2112: been supplied explicitly by -pc_fieldsplit_%d_fields
2113: . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields)
2114: . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur> - type of relaxation or factorization splitting
2115: . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11
2116: . -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero or negative diagonal and uses Schur complement with no preconditioner as the solver
2118: - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_
2119: for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields
2121: Notes:
2122: Use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS()
2123: to define a field by an arbitrary collection of entries.
2125: If no fields are set the default is used. The fields are defined by entries strided by bs,
2126: beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(),
2127: if this is not called the block size defaults to the blocksize of the second matrix passed
2128: to KSPSetOperators()/PCSetOperators().
2130: $ For the Schur complement preconditioner if J = ( A00 A01 )
2131: $ ( A10 A11 )
2132: $ the preconditioner using full factorization is
2133: $ ( I -ksp(A00) A01 ) ( inv(A00) 0 ) ( I 0 )
2134: $ ( 0 I ) ( 0 ksp(S) ) ( -A10 ksp(A00) I )
2135: where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. S is the Schur complement
2136: $ S = A11 - A10 ksp(A00) A01
2137: 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
2138: in providing the SECOND split or 1 if not give). For PCFieldSplitGetKSP() when field number is 0,
2139: it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default
2140: A11 is used to construct a preconditioner for S, use PCFieldSplitSetSchurPre() to turn on or off this
2141: option. You can use the preconditioner PCLSC to precondition the Schur complement with -fieldsplit_1_pc_type lsc.
2142: When option -fieldsplit_schur_precondition selfp is given, an approximation to S is assembled --
2143: Sp = A11 - A10 inv(diag(A00)) A01, which has type AIJ and can be used with a variety of preconditioners
2144: (e.g., -fieldsplit_1_pc_type asm).
2145: The factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above,
2146: diag gives
2147: $ ( inv(A00) 0 )
2148: $ ( 0 -ksp(S) )
2149: 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
2150: $ ( A00 0 )
2151: $ ( A10 S )
2152: where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of
2153: $ ( A00 A01 )
2154: $ ( 0 S )
2155: where again the inverses of A00 and S are applied using KSPs.
2157: If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS
2158: is used automatically for a second block.
2160: The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1.
2161: Generally it should be used with the AIJ format.
2163: The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see,
2164: for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT
2165: inside a smoother resulting in "Distributive Smoothers".
2167: Concepts: physics based preconditioners, block preconditioners
2169: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC,
2170: PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre()
2171: M*/
2175: PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc)
2176: {
2178: PC_FieldSplit *jac;
2181: PetscNewLog(pc,&jac);
2183: jac->bs = -1;
2184: jac->nsplits = 0;
2185: jac->type = PC_COMPOSITE_MULTIPLICATIVE;
2186: jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */
2187: jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL;
2188: jac->dm_splits = PETSC_TRUE;
2190: pc->data = (void*)jac;
2192: pc->ops->apply = PCApply_FieldSplit;
2193: pc->ops->applytranspose = PCApplyTranspose_FieldSplit;
2194: pc->ops->setup = PCSetUp_FieldSplit;
2195: pc->ops->reset = PCReset_FieldSplit;
2196: pc->ops->destroy = PCDestroy_FieldSplit;
2197: pc->ops->setfromoptions = PCSetFromOptions_FieldSplit;
2198: pc->ops->view = PCView_FieldSplit;
2199: pc->ops->applyrichardson = 0;
2201: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2202: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);
2203: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);
2204: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);
2205: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);
2206: return(0);
2207: }