Actual source code: fieldsplit.c
petsc-3.11.4 2019-09-28
3: #include <petsc/private/pcimpl.h>
4: #include <petsc/private/kspimpl.h> /* This is needed to provide the appropriate PETSC_EXTERN for KSP_Solve_FS ....*/
5: #include <petscdm.h>
7: const char *const PCFieldSplitSchurPreTypes[] = {"SELF","SELFP","A11","USER","FULL","PCFieldSplitSchurPreType","PC_FIELDSPLIT_SCHUR_PRE_",0};
8: const char *const PCFieldSplitSchurFactTypes[] = {"DIAG","LOWER","UPPER","FULL","PCFieldSplitSchurFactType","PC_FIELDSPLIT_SCHUR_FACT_",0};
10: PetscLogEvent KSP_Solve_FS_0,KSP_Solve_FS_1,KSP_Solve_FS_S,KSP_Solve_FS_U,KSP_Solve_FS_L,KSP_Solve_FS_2,KSP_Solve_FS_3,KSP_Solve_FS_4;
12: typedef struct _PC_FieldSplitLink *PC_FieldSplitLink;
13: struct _PC_FieldSplitLink {
14: KSP ksp;
15: Vec x,y,z;
16: char *splitname;
17: PetscInt nfields;
18: PetscInt *fields,*fields_col;
19: VecScatter sctx;
20: IS is,is_col;
21: PC_FieldSplitLink next,previous;
22: PetscLogEvent event;
23: };
25: typedef struct {
26: PCCompositeType type;
27: PetscBool defaultsplit; /* Flag for a system with a set of 'k' scalar fields with the same layout (and bs = k) */
28: PetscBool splitdefined; /* Flag is set after the splits have been defined, to prevent more splits from being added */
29: PetscInt bs; /* Block size for IS and Mat structures */
30: PetscInt nsplits; /* Number of field divisions defined */
31: Vec *x,*y,w1,w2;
32: Mat *mat; /* The diagonal block for each split */
33: Mat *pmat; /* The preconditioning diagonal block for each split */
34: Mat *Afield; /* The rows of the matrix associated with each split */
35: PetscBool issetup;
37: /* Only used when Schur complement preconditioning is used */
38: Mat B; /* The (0,1) block */
39: Mat C; /* The (1,0) block */
40: Mat schur; /* The Schur complement S = A11 - A10 A00^{-1} A01, the KSP here, kspinner, is H_1 in [El08] */
41: Mat schurp; /* Assembled approximation to S built by MatSchurComplement to be used as a preconditioning matrix when solving with S */
42: Mat schur_user; /* User-provided preconditioning matrix for the Schur complement */
43: PCFieldSplitSchurPreType schurpre; /* Determines which preconditioning matrix is used for the Schur complement */
44: PCFieldSplitSchurFactType schurfactorization;
45: KSP kspschur; /* The solver for S */
46: KSP kspupper; /* The solver for A in the upper diagonal part of the factorization (H_2 in [El08]) */
47: PetscScalar schurscale; /* Scaling factor for the Schur complement solution with DIAG factorization */
49: /* Only used when Golub-Kahan bidiagonalization preconditioning is used */
50: Mat H; /* The modified matrix H = A00 + nu*A01*A01' */
51: PetscReal gkbtol; /* Stopping tolerance for lower bound estimate */
52: PetscInt gkbdelay; /* The delay window for the stopping criterion */
53: PetscReal gkbnu; /* Parameter for augmented Lagrangian H = A + nu*A01*A01' */
54: PetscInt gkbmaxit; /* Maximum number of iterations for outer loop */
55: PetscBool gkbmonitor; /* Monitor for gkb iterations and the lower bound error */
56: PetscViewer gkbviewer; /* Viewer context for gkbmonitor */
57: Vec u,v,d,Hu; /* Work vectors for the GKB algorithm */
58: PetscScalar *vecz; /* Contains intermediate values, eg for lower bound */
60: PC_FieldSplitLink head;
61: PetscBool isrestrict; /* indicates PCFieldSplitRestrictIS() has been last called on this object, hack */
62: PetscBool suboptionsset; /* Indicates that the KSPSetFromOptions() has been called on the sub-KSPs */
63: PetscBool dm_splits; /* Whether to use DMCreateFieldDecomposition() whenever possible */
64: PetscBool diag_use_amat; /* Whether to extract diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */
65: PetscBool offdiag_use_amat; /* Whether to extract off-diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */
66: PetscBool detect; /* Whether to form 2-way split by finding zero diagonal entries */
67: } PC_FieldSplit;
69: /*
70: Notes:
71: there is no particular reason that pmat, x, and y are stored as arrays in PC_FieldSplit instead of
72: inside PC_FieldSplitLink, just historical. If you want to be able to add new fields after already using the
73: PC you could change this.
74: */
76: /* This helper is so that setting a user-provided preconditioning matrix is orthogonal to choosing to use it. This way the
77: * application-provided FormJacobian can provide this matrix without interfering with the user's (command-line) choices. */
78: static Mat FieldSplitSchurPre(PC_FieldSplit *jac)
79: {
80: switch (jac->schurpre) {
81: case PC_FIELDSPLIT_SCHUR_PRE_SELF: return jac->schur;
82: case PC_FIELDSPLIT_SCHUR_PRE_SELFP: return jac->schurp;
83: case PC_FIELDSPLIT_SCHUR_PRE_A11: return jac->pmat[1];
84: case PC_FIELDSPLIT_SCHUR_PRE_FULL: /* We calculate this and store it in schur_user */
85: case PC_FIELDSPLIT_SCHUR_PRE_USER: /* Use a user-provided matrix if it is given, otherwise diagonal block */
86: default:
87: return jac->schur_user ? jac->schur_user : jac->pmat[1];
88: }
89: }
92: #include <petscdraw.h>
93: static PetscErrorCode PCView_FieldSplit(PC pc,PetscViewer viewer)
94: {
95: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
96: PetscErrorCode ierr;
97: PetscBool iascii,isdraw;
98: PetscInt i,j;
99: PC_FieldSplitLink ilink = jac->head;
102: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
103: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
104: if (iascii) {
105: if (jac->bs > 0) {
106: PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs);
107: } else {
108: PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D\n",PCCompositeTypes[jac->type],jac->nsplits);
109: }
110: if (pc->useAmat) {
111: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");
112: }
113: if (jac->diag_use_amat) {
114: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for diagonal blocks\n");
115: }
116: if (jac->offdiag_use_amat) {
117: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for off-diagonal blocks\n");
118: }
119: PetscViewerASCIIPrintf(viewer," Solver info for each split is in the following KSP objects:\n");
120: for (i=0; i<jac->nsplits; i++) {
121: if (ilink->fields) {
122: PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);
123: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
124: for (j=0; j<ilink->nfields; j++) {
125: if (j > 0) {
126: PetscViewerASCIIPrintf(viewer,",");
127: }
128: PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);
129: }
130: PetscViewerASCIIPrintf(viewer,"\n");
131: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
132: } else {
133: PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);
134: }
135: KSPView(ilink->ksp,viewer);
136: ilink = ilink->next;
137: }
138: }
140: if (isdraw) {
141: PetscDraw draw;
142: PetscReal x,y,w,wd;
144: PetscViewerDrawGetDraw(viewer,0,&draw);
145: PetscDrawGetCurrentPoint(draw,&x,&y);
146: w = 2*PetscMin(1.0 - x,x);
147: wd = w/(jac->nsplits + 1);
148: x = x - wd*(jac->nsplits-1)/2.0;
149: for (i=0; i<jac->nsplits; i++) {
150: PetscDrawPushCurrentPoint(draw,x,y);
151: KSPView(ilink->ksp,viewer);
152: PetscDrawPopCurrentPoint(draw);
153: x += wd;
154: ilink = ilink->next;
155: }
156: }
157: return(0);
158: }
160: static PetscErrorCode PCView_FieldSplit_Schur(PC pc,PetscViewer viewer)
161: {
162: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
163: PetscErrorCode ierr;
164: PetscBool iascii,isdraw;
165: PetscInt i,j;
166: PC_FieldSplitLink ilink = jac->head;
167: MatSchurComplementAinvType atype;
170: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
171: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
172: if (iascii) {
173: if (jac->bs > 0) {
174: PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, blocksize = %D, factorization %s\n",jac->bs,PCFieldSplitSchurFactTypes[jac->schurfactorization]);
175: } else {
176: PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, factorization %s\n",PCFieldSplitSchurFactTypes[jac->schurfactorization]);
177: }
178: if (pc->useAmat) {
179: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");
180: }
181: switch (jac->schurpre) {
182: case PC_FIELDSPLIT_SCHUR_PRE_SELF:
183: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from S itself\n");
184: break;
185: case PC_FIELDSPLIT_SCHUR_PRE_SELFP:
186: MatSchurComplementGetAinvType(jac->schur,&atype);
187: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from Sp, an assembled approximation to S, which uses A00's %sdiagonal's inverse\n",atype == MAT_SCHUR_COMPLEMENT_AINV_DIAG ? "" : (atype == MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG ? "block " : "lumped "));break;
188: case PC_FIELDSPLIT_SCHUR_PRE_A11:
189: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n");
190: break;
191: case PC_FIELDSPLIT_SCHUR_PRE_FULL:
192: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from the exact Schur complement\n");
193: break;
194: case PC_FIELDSPLIT_SCHUR_PRE_USER:
195: if (jac->schur_user) {
196: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from user provided matrix\n");
197: } else {
198: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n");
199: }
200: break;
201: default:
202: SETERRQ1(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Invalid Schur preconditioning type: %d", jac->schurpre);
203: }
204: PetscViewerASCIIPrintf(viewer," Split info:\n");
205: PetscViewerASCIIPushTab(viewer);
206: for (i=0; i<jac->nsplits; i++) {
207: if (ilink->fields) {
208: PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);
209: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
210: for (j=0; j<ilink->nfields; j++) {
211: if (j > 0) {
212: PetscViewerASCIIPrintf(viewer,",");
213: }
214: PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);
215: }
216: PetscViewerASCIIPrintf(viewer,"\n");
217: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
218: } else {
219: PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);
220: }
221: ilink = ilink->next;
222: }
223: PetscViewerASCIIPrintf(viewer,"KSP solver for A00 block\n");
224: PetscViewerASCIIPushTab(viewer);
225: if (jac->head) {
226: KSPView(jac->head->ksp,viewer);
227: } else {PetscViewerASCIIPrintf(viewer," not yet available\n");}
228: PetscViewerASCIIPopTab(viewer);
229: if (jac->head && jac->kspupper != jac->head->ksp) {
230: PetscViewerASCIIPrintf(viewer,"KSP solver for upper A00 in upper triangular factor \n");
231: PetscViewerASCIIPushTab(viewer);
232: if (jac->kspupper) {KSPView(jac->kspupper,viewer);}
233: else {PetscViewerASCIIPrintf(viewer," not yet available\n");}
234: PetscViewerASCIIPopTab(viewer);
235: }
236: PetscViewerASCIIPrintf(viewer,"KSP solver for S = A11 - A10 inv(A00) A01 \n");
237: PetscViewerASCIIPushTab(viewer);
238: if (jac->kspschur) {
239: KSPView(jac->kspschur,viewer);
240: } else {
241: PetscViewerASCIIPrintf(viewer," not yet available\n");
242: }
243: PetscViewerASCIIPopTab(viewer);
244: PetscViewerASCIIPopTab(viewer);
245: } else if (isdraw && jac->head) {
246: PetscDraw draw;
247: PetscReal x,y,w,wd,h;
248: PetscInt cnt = 2;
249: char str[32];
251: PetscViewerDrawGetDraw(viewer,0,&draw);
252: PetscDrawGetCurrentPoint(draw,&x,&y);
253: if (jac->kspupper != jac->head->ksp) cnt++;
254: w = 2*PetscMin(1.0 - x,x);
255: wd = w/(cnt + 1);
257: PetscSNPrintf(str,32,"Schur fact. %s",PCFieldSplitSchurFactTypes[jac->schurfactorization]);
258: PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
259: y -= h;
260: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_USER && !jac->schur_user) {
261: PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[PC_FIELDSPLIT_SCHUR_PRE_A11]);
262: } else {
263: PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[jac->schurpre]);
264: }
265: PetscDrawStringBoxed(draw,x+wd*(cnt-1)/2.0,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
266: y -= h;
267: x = x - wd*(cnt-1)/2.0;
269: PetscDrawPushCurrentPoint(draw,x,y);
270: KSPView(jac->head->ksp,viewer);
271: PetscDrawPopCurrentPoint(draw);
272: if (jac->kspupper != jac->head->ksp) {
273: x += wd;
274: PetscDrawPushCurrentPoint(draw,x,y);
275: KSPView(jac->kspupper,viewer);
276: PetscDrawPopCurrentPoint(draw);
277: }
278: x += wd;
279: PetscDrawPushCurrentPoint(draw,x,y);
280: KSPView(jac->kspschur,viewer);
281: PetscDrawPopCurrentPoint(draw);
282: }
283: return(0);
284: }
286: static PetscErrorCode PCView_FieldSplit_GKB(PC pc,PetscViewer viewer)
287: {
288: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
289: PetscErrorCode ierr;
290: PetscBool iascii,isdraw;
291: PetscInt i,j;
292: PC_FieldSplitLink ilink = jac->head;
295: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
296: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
297: if (iascii) {
298: if (jac->bs > 0) {
299: PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs);
300: } else {
301: PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D\n",PCCompositeTypes[jac->type],jac->nsplits);
302: }
303: if (pc->useAmat) {
304: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");
305: }
306: if (jac->diag_use_amat) {
307: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for diagonal blocks\n");
308: }
309: if (jac->offdiag_use_amat) {
310: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for off-diagonal blocks\n");
311: }
313: PetscViewerASCIIPrintf(viewer," Stopping tolerance=%.1e, delay in error estimate=%D, maximum iterations=%D\n",jac->gkbtol,jac->gkbdelay,jac->gkbmaxit);
314: PetscViewerASCIIPrintf(viewer," Solver info for H = A00 + nu*A01*A01' matrix:\n");
315: PetscViewerASCIIPushTab(viewer);
317: if (ilink->fields) {
318: PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",0);
319: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
320: for (j=0; j<ilink->nfields; j++) {
321: if (j > 0) {
322: PetscViewerASCIIPrintf(viewer,",");
323: }
324: PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);
325: }
326: PetscViewerASCIIPrintf(viewer,"\n");
327: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
328: } else {
329: PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",0);
330: }
331: KSPView(ilink->ksp,viewer);
333: PetscViewerASCIIPopTab(viewer);
334: }
336: if (isdraw) {
337: PetscDraw draw;
338: PetscReal x,y,w,wd;
340: PetscViewerDrawGetDraw(viewer,0,&draw);
341: PetscDrawGetCurrentPoint(draw,&x,&y);
342: w = 2*PetscMin(1.0 - x,x);
343: wd = w/(jac->nsplits + 1);
344: x = x - wd*(jac->nsplits-1)/2.0;
345: for (i=0; i<jac->nsplits; i++) {
346: PetscDrawPushCurrentPoint(draw,x,y);
347: KSPView(ilink->ksp,viewer);
348: PetscDrawPopCurrentPoint(draw);
349: x += wd;
350: ilink = ilink->next;
351: }
352: }
353: return(0);
354: }
357: /* Precondition: jac->bs is set to a meaningful value */
358: static PetscErrorCode PCFieldSplitSetRuntimeSplits_Private(PC pc)
359: {
361: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
362: PetscInt i,nfields,*ifields,nfields_col,*ifields_col;
363: PetscBool flg,flg_col;
364: char optionname[128],splitname[8],optionname_col[128];
367: PetscMalloc1(jac->bs,&ifields);
368: PetscMalloc1(jac->bs,&ifields_col);
369: for (i=0,flg=PETSC_TRUE;; i++) {
370: PetscSNPrintf(splitname,sizeof(splitname),"%D",i);
371: PetscSNPrintf(optionname,sizeof(optionname),"-pc_fieldsplit_%D_fields",i);
372: PetscSNPrintf(optionname_col,sizeof(optionname_col),"-pc_fieldsplit_%D_fields_col",i);
373: nfields = jac->bs;
374: nfields_col = jac->bs;
375: PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,optionname,ifields,&nfields,&flg);
376: PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,optionname_col,ifields_col,&nfields_col,&flg_col);
377: if (!flg) break;
378: else if (flg && !flg_col) {
379: if (!nfields) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields");
380: PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields);
381: } else {
382: if (!nfields || !nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields");
383: if (nfields != nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Number of row and column fields must match");
384: PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields_col);
385: }
386: }
387: if (i > 0) {
388: /* Makes command-line setting of splits take precedence over setting them in code.
389: Otherwise subsequent calls to PCFieldSplitSetIS() or PCFieldSplitSetFields() would
390: create new splits, which would probably not be what the user wanted. */
391: jac->splitdefined = PETSC_TRUE;
392: }
393: PetscFree(ifields);
394: PetscFree(ifields_col);
395: return(0);
396: }
398: static PetscErrorCode PCFieldSplitSetDefaults(PC pc)
399: {
400: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
401: PetscErrorCode ierr;
402: PC_FieldSplitLink ilink = jac->head;
403: PetscBool fieldsplit_default = PETSC_FALSE,coupling = PETSC_FALSE;
404: PetscInt i;
407: /*
408: Kinda messy, but at least this now uses DMCreateFieldDecomposition().
409: Should probably be rewritten.
410: */
411: if (!ilink) {
412: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_coupling",&coupling,NULL);
413: if (pc->dm && jac->dm_splits && !jac->detect && !coupling) {
414: PetscInt numFields, f, i, j;
415: char **fieldNames;
416: IS *fields;
417: DM *dms;
418: DM subdm[128];
419: PetscBool flg;
421: DMCreateFieldDecomposition(pc->dm, &numFields, &fieldNames, &fields, &dms);
422: /* Allow the user to prescribe the splits */
423: for (i = 0, flg = PETSC_TRUE;; i++) {
424: PetscInt ifields[128];
425: IS compField;
426: char optionname[128], splitname[8];
427: PetscInt nfields = numFields;
429: PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%D_fields", i);
430: PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix, optionname, ifields, &nfields, &flg);
431: if (!flg) break;
432: if (numFields > 128) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot currently support %d > 128 fields", numFields);
433: DMCreateSubDM(pc->dm, nfields, ifields, &compField, &subdm[i]);
434: if (nfields == 1) {
435: PCFieldSplitSetIS(pc, fieldNames[ifields[0]], compField);
436: } else {
437: PetscSNPrintf(splitname, sizeof(splitname), "%D", i);
438: PCFieldSplitSetIS(pc, splitname, compField);
439: }
440: ISDestroy(&compField);
441: for (j = 0; j < nfields; ++j) {
442: f = ifields[j];
443: PetscFree(fieldNames[f]);
444: ISDestroy(&fields[f]);
445: }
446: }
447: if (i == 0) {
448: for (f = 0; f < numFields; ++f) {
449: PCFieldSplitSetIS(pc, fieldNames[f], fields[f]);
450: PetscFree(fieldNames[f]);
451: ISDestroy(&fields[f]);
452: }
453: } else {
454: for (j=0; j<numFields; j++) {
455: DMDestroy(dms+j);
456: }
457: PetscFree(dms);
458: PetscMalloc1(i, &dms);
459: for (j = 0; j < i; ++j) dms[j] = subdm[j];
460: }
461: PetscFree(fieldNames);
462: PetscFree(fields);
463: if (dms) {
464: PetscInfo(pc, "Setting up physics based fieldsplit preconditioner using the embedded DM\n");
465: for (ilink = jac->head, i = 0; ilink; ilink = ilink->next, ++i) {
466: const char *prefix;
467: PetscObjectGetOptionsPrefix((PetscObject)(ilink->ksp),&prefix);
468: PetscObjectSetOptionsPrefix((PetscObject)(dms[i]), prefix);
469: KSPSetDM(ilink->ksp, dms[i]);
470: KSPSetDMActive(ilink->ksp, PETSC_FALSE);
471: {
472: PetscErrorCode (*func)(KSP,Mat,Mat,void*);
473: void *ctx;
475: DMKSPGetComputeOperators(pc->dm, &func, &ctx);
476: DMKSPSetComputeOperators(dms[i], func, ctx);
477: }
478: PetscObjectIncrementTabLevel((PetscObject)dms[i],(PetscObject)ilink->ksp,0);
479: DMDestroy(&dms[i]);
480: }
481: PetscFree(dms);
482: }
483: } else {
484: if (jac->bs <= 0) {
485: if (pc->pmat) {
486: MatGetBlockSize(pc->pmat,&jac->bs);
487: } else jac->bs = 1;
488: }
490: if (jac->detect) {
491: IS zerodiags,rest;
492: PetscInt nmin,nmax;
494: MatGetOwnershipRange(pc->mat,&nmin,&nmax);
495: MatFindZeroDiagonals(pc->mat,&zerodiags);
496: ISComplement(zerodiags,nmin,nmax,&rest);
497: PCFieldSplitSetIS(pc,"0",rest);
498: PCFieldSplitSetIS(pc,"1",zerodiags);
499: ISDestroy(&zerodiags);
500: ISDestroy(&rest);
501: } else if (coupling) {
502: IS coupling,rest;
503: PetscInt nmin,nmax;
505: MatGetOwnershipRange(pc->mat,&nmin,&nmax);
506: MatFindOffBlockDiagonalEntries(pc->mat,&coupling);
507: ISCreateStride(PetscObjectComm((PetscObject)pc->mat),nmax-nmin,nmin,1,&rest);
508: ISSetIdentity(rest);
509: PCFieldSplitSetIS(pc,"0",rest);
510: PCFieldSplitSetIS(pc,"1",coupling);
511: ISDestroy(&coupling);
512: ISDestroy(&rest);
513: } else {
514: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&fieldsplit_default,NULL);
515: if (!fieldsplit_default) {
516: /* Allow user to set fields from command line, if bs was known at the time of PCSetFromOptions_FieldSplit()
517: then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */
518: PCFieldSplitSetRuntimeSplits_Private(pc);
519: if (jac->splitdefined) {PetscInfo(pc,"Splits defined using the options database\n");}
520: }
521: if ((fieldsplit_default || !jac->splitdefined) && !jac->isrestrict) {
522: PetscInfo(pc,"Using default splitting of fields\n");
523: for (i=0; i<jac->bs; i++) {
524: char splitname[8];
525: PetscSNPrintf(splitname,sizeof(splitname),"%D",i);
526: PCFieldSplitSetFields(pc,splitname,1,&i,&i);
527: }
528: jac->defaultsplit = PETSC_TRUE;
529: }
530: }
531: }
532: } else if (jac->nsplits == 1) {
533: if (ilink->is) {
534: IS is2;
535: PetscInt nmin,nmax;
537: MatGetOwnershipRange(pc->mat,&nmin,&nmax);
538: ISComplement(ilink->is,nmin,nmax,&is2);
539: PCFieldSplitSetIS(pc,"1",is2);
540: ISDestroy(&is2);
541: } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Must provide at least two sets of fields to PCFieldSplit()");
542: }
544: if (jac->nsplits < 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Unhandled case, must have at least two fields, not %d", jac->nsplits);
545: return(0);
546: }
548: static PetscErrorCode MatGolubKahanComputeExplicitOperator(Mat A,Mat B,Mat C,Mat *H,PetscReal gkbnu)
549: {
550: PetscErrorCode ierr;
551: Mat BT,T;
552: PetscReal nrmT,nrmB;
555: MatHermitianTranspose(C,MAT_INITIAL_MATRIX,&T); /* Test if augmented matrix is symmetric */
556: MatAXPY(T,-1.0,B,DIFFERENT_NONZERO_PATTERN);
557: MatNorm(T,NORM_1,&nrmT);
558: MatNorm(B,NORM_1,&nrmB);
559: if (nrmB > 0) {
560: if (nrmT/nrmB >= PETSC_SMALL) {
561: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix is not symmetric/hermitian, GKB is not applicable.");
562: }
563: }
564: /* Compute augmented Lagrangian matrix H = A00 + nu*A01*A01'. This corresponds to */
565: /* setting N := 1/nu*I in [Ar13]. */
566: MatHermitianTranspose(B,MAT_INITIAL_MATRIX,&BT);
567: MatMatMult(B,BT,MAT_INITIAL_MATRIX,PETSC_DEFAULT,H); /* H = A01*A01' */
568: MatAYPX(*H,gkbnu,A,DIFFERENT_NONZERO_PATTERN); /* H = A00 + nu*A01*A01' */
570: MatDestroy(&BT);
571: MatDestroy(&T);
572: return(0);
573: }
575: PETSC_EXTERN PetscErrorCode PetscOptionsFindPairPrefix_Private(PetscOptions,const char pre[], const char name[],const char *value[],PetscBool *flg);
577: static PetscErrorCode PCSetUp_FieldSplit(PC pc)
578: {
579: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
580: PetscErrorCode ierr;
581: PC_FieldSplitLink ilink;
582: PetscInt i,nsplit;
583: PetscBool sorted, sorted_col;
586: pc->failedreason = PC_NOERROR;
587: PCFieldSplitSetDefaults(pc);
588: nsplit = jac->nsplits;
589: ilink = jac->head;
591: /* get the matrices for each split */
592: if (!jac->issetup) {
593: PetscInt rstart,rend,nslots,bs;
595: jac->issetup = PETSC_TRUE;
597: /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */
598: if (jac->defaultsplit || !ilink->is) {
599: if (jac->bs <= 0) jac->bs = nsplit;
600: }
601: bs = jac->bs;
602: MatGetOwnershipRange(pc->pmat,&rstart,&rend);
603: nslots = (rend - rstart)/bs;
604: for (i=0; i<nsplit; i++) {
605: if (jac->defaultsplit) {
606: ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+i,nsplit,&ilink->is);
607: ISDuplicate(ilink->is,&ilink->is_col);
608: } else if (!ilink->is) {
609: if (ilink->nfields > 1) {
610: PetscInt *ii,*jj,j,k,nfields = ilink->nfields,*fields = ilink->fields,*fields_col = ilink->fields_col;
611: PetscMalloc1(ilink->nfields*nslots,&ii);
612: PetscMalloc1(ilink->nfields*nslots,&jj);
613: for (j=0; j<nslots; j++) {
614: for (k=0; k<nfields; k++) {
615: ii[nfields*j + k] = rstart + bs*j + fields[k];
616: jj[nfields*j + k] = rstart + bs*j + fields_col[k];
617: }
618: }
619: ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,ii,PETSC_OWN_POINTER,&ilink->is);
620: ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,jj,PETSC_OWN_POINTER,&ilink->is_col);
621: ISSetBlockSize(ilink->is, nfields);
622: ISSetBlockSize(ilink->is_col, nfields);
623: } else {
624: ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields[0],bs,&ilink->is);
625: ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields_col[0],bs,&ilink->is_col);
626: }
627: }
628: ISSorted(ilink->is,&sorted);
629: if (ilink->is_col) { ISSorted(ilink->is_col,&sorted_col); }
630: if (!sorted || !sorted_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Fields must be sorted when creating split");
631: ilink = ilink->next;
632: }
633: }
635: ilink = jac->head;
636: if (!jac->pmat) {
637: Vec xtmp;
639: MatCreateVecs(pc->pmat,&xtmp,NULL);
640: PetscMalloc1(nsplit,&jac->pmat);
641: PetscMalloc2(nsplit,&jac->x,nsplit,&jac->y);
642: for (i=0; i<nsplit; i++) {
643: MatNullSpace sp;
645: /* Check for preconditioning matrix attached to IS */
646: PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &jac->pmat[i]);
647: if (jac->pmat[i]) {
648: PetscObjectReference((PetscObject) jac->pmat[i]);
649: if (jac->type == PC_COMPOSITE_SCHUR) {
650: jac->schur_user = jac->pmat[i];
652: PetscObjectReference((PetscObject) jac->schur_user);
653: }
654: } else {
655: const char *prefix;
656: MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->pmat[i]);
657: KSPGetOptionsPrefix(ilink->ksp,&prefix);
658: MatSetOptionsPrefix(jac->pmat[i],prefix);
659: MatViewFromOptions(jac->pmat[i],NULL,"-mat_view");
660: }
661: /* create work vectors for each split */
662: MatCreateVecs(jac->pmat[i],&jac->x[i],&jac->y[i]);
663: ilink->x = jac->x[i]; ilink->y = jac->y[i]; ilink->z = NULL;
664: /* compute scatter contexts needed by multiplicative versions and non-default splits */
665: VecScatterCreate(xtmp,ilink->is,jac->x[i],NULL,&ilink->sctx);
666: PetscObjectQuery((PetscObject) ilink->is, "nearnullspace", (PetscObject*) &sp);
667: if (sp) {
668: MatSetNearNullSpace(jac->pmat[i], sp);
669: }
670: ilink = ilink->next;
671: }
672: VecDestroy(&xtmp);
673: } else {
674: MatReuse scall;
675: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
676: for (i=0; i<nsplit; i++) {
677: MatDestroy(&jac->pmat[i]);
678: }
679: scall = MAT_INITIAL_MATRIX;
680: } else scall = MAT_REUSE_MATRIX;
682: for (i=0; i<nsplit; i++) {
683: Mat pmat;
685: /* Check for preconditioning matrix attached to IS */
686: PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &pmat);
687: if (!pmat) {
688: MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is_col,scall,&jac->pmat[i]);
689: }
690: ilink = ilink->next;
691: }
692: }
693: if (jac->diag_use_amat) {
694: ilink = jac->head;
695: if (!jac->mat) {
696: PetscMalloc1(nsplit,&jac->mat);
697: for (i=0; i<nsplit; i++) {
698: MatCreateSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->mat[i]);
699: ilink = ilink->next;
700: }
701: } else {
702: MatReuse scall;
703: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
704: for (i=0; i<nsplit; i++) {
705: MatDestroy(&jac->mat[i]);
706: }
707: scall = MAT_INITIAL_MATRIX;
708: } else scall = MAT_REUSE_MATRIX;
710: for (i=0; i<nsplit; i++) {
711: if (jac->mat[i]) {MatCreateSubMatrix(pc->mat,ilink->is,ilink->is_col,scall,&jac->mat[i]);}
712: ilink = ilink->next;
713: }
714: }
715: } else {
716: jac->mat = jac->pmat;
717: }
719: /* Check for null space attached to IS */
720: ilink = jac->head;
721: for (i=0; i<nsplit; i++) {
722: MatNullSpace sp;
724: PetscObjectQuery((PetscObject) ilink->is, "nullspace", (PetscObject*) &sp);
725: if (sp) {
726: MatSetNullSpace(jac->mat[i], sp);
727: }
728: ilink = ilink->next;
729: }
731: if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR && jac->type != PC_COMPOSITE_GKB) {
732: /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */
733: /* FIXME: Can/should we reuse jac->mat whenever (jac->diag_use_amat) is true? */
734: ilink = jac->head;
735: if (nsplit == 2 && jac->type == PC_COMPOSITE_MULTIPLICATIVE) {
736: /* 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 */
737: if (!jac->Afield) {
738: PetscCalloc1(nsplit,&jac->Afield);
739: if (jac->offdiag_use_amat) {
740: MatCreateSubMatrix(pc->mat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);
741: } else {
742: MatCreateSubMatrix(pc->pmat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);
743: }
744: } else {
745: MatReuse scall;
746: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
747: for (i=0; i<nsplit; i++) {
748: MatDestroy(&jac->Afield[1]);
749: }
750: scall = MAT_INITIAL_MATRIX;
751: } else scall = MAT_REUSE_MATRIX;
753: if (jac->offdiag_use_amat) {
754: MatCreateSubMatrix(pc->mat,ilink->next->is,ilink->is,scall,&jac->Afield[1]);
755: } else {
756: MatCreateSubMatrix(pc->pmat,ilink->next->is,ilink->is,scall,&jac->Afield[1]);
757: }
758: }
759: } else {
760: if (!jac->Afield) {
761: PetscMalloc1(nsplit,&jac->Afield);
762: for (i=0; i<nsplit; i++) {
763: if (jac->offdiag_use_amat) {
764: MatCreateSubMatrix(pc->mat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);
765: } else {
766: MatCreateSubMatrix(pc->pmat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);
767: }
768: ilink = ilink->next;
769: }
770: } else {
771: MatReuse scall;
772: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
773: for (i=0; i<nsplit; i++) {
774: MatDestroy(&jac->Afield[i]);
775: }
776: scall = MAT_INITIAL_MATRIX;
777: } else scall = MAT_REUSE_MATRIX;
779: for (i=0; i<nsplit; i++) {
780: if (jac->offdiag_use_amat) {
781: MatCreateSubMatrix(pc->mat,ilink->is,NULL,scall,&jac->Afield[i]);
782: } else {
783: MatCreateSubMatrix(pc->pmat,ilink->is,NULL,scall,&jac->Afield[i]);
784: }
785: ilink = ilink->next;
786: }
787: }
788: }
789: }
791: if (jac->type == PC_COMPOSITE_SCHUR) {
792: IS ccis;
793: PetscBool isspd;
794: PetscInt rstart,rend;
795: char lscname[256];
796: PetscObject LSC_L;
798: if (nsplit != 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"To use Schur complement preconditioner you must have exactly 2 fields");
800: /* If pc->mat is SPD, don't scale by -1 the Schur complement */
801: if (jac->schurscale == (PetscScalar)-1.0) {
802: MatGetOption(pc->pmat,MAT_SPD,&isspd);
803: jac->schurscale = (isspd == PETSC_TRUE) ? 1.0 : -1.0;
804: }
806: /* When extracting off-diagonal submatrices, we take complements from this range */
807: MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);
809: /* need to handle case when one is resetting up the preconditioner */
810: if (jac->schur) {
811: KSP kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper;
813: MatSchurComplementGetKSP(jac->schur, &kspInner);
814: ilink = jac->head;
815: ISComplement(ilink->is_col,rstart,rend,&ccis);
816: if (jac->offdiag_use_amat) {
817: MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);
818: } else {
819: MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);
820: }
821: ISDestroy(&ccis);
822: ilink = ilink->next;
823: ISComplement(ilink->is_col,rstart,rend,&ccis);
824: if (jac->offdiag_use_amat) {
825: MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);
826: } else {
827: MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);
828: }
829: ISDestroy(&ccis);
830: MatSchurComplementUpdateSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);
831: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
832: MatDestroy(&jac->schurp);
833: MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);
834: }
835: if (kspA != kspInner) {
836: KSPSetOperators(kspA,jac->mat[0],jac->pmat[0]);
837: }
838: if (kspUpper != kspA) {
839: KSPSetOperators(kspUpper,jac->mat[0],jac->pmat[0]);
840: }
841: KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));
842: } else {
843: const char *Dprefix;
844: char schurprefix[256], schurmatprefix[256];
845: char schurtestoption[256];
846: MatNullSpace sp;
847: PetscBool flg;
848: KSP kspt;
850: /* extract the A01 and A10 matrices */
851: ilink = jac->head;
852: ISComplement(ilink->is_col,rstart,rend,&ccis);
853: if (jac->offdiag_use_amat) {
854: MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
855: } else {
856: MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
857: }
858: ISDestroy(&ccis);
859: ilink = ilink->next;
860: ISComplement(ilink->is_col,rstart,rend,&ccis);
861: if (jac->offdiag_use_amat) {
862: MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
863: } else {
864: MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
865: }
866: ISDestroy(&ccis);
868: /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */
869: MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur);
870: MatSetType(jac->schur,MATSCHURCOMPLEMENT);
871: MatSchurComplementSetSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);
872: PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
873: MatSetOptionsPrefix(jac->schur,schurmatprefix);
874: MatSchurComplementGetKSP(jac->schur,&kspt);
875: KSPSetOptionsPrefix(kspt,schurmatprefix);
877: /* Note: this is not true in general */
878: MatGetNullSpace(jac->mat[1], &sp);
879: if (sp) {
880: MatSetNullSpace(jac->schur, sp);
881: }
883: PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);
884: PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
885: if (flg) {
886: DM dmInner;
887: KSP kspInner;
888: PC pcInner;
890: MatSchurComplementGetKSP(jac->schur, &kspInner);
891: KSPReset(kspInner);
892: KSPSetOperators(kspInner,jac->mat[0],jac->pmat[0]);
893: PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
894: /* Indent this deeper to emphasize the "inner" nature of this solver. */
895: PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject) pc, 2);
896: PetscObjectIncrementTabLevel((PetscObject)kspInner->pc, (PetscObject) pc, 2);
897: KSPSetOptionsPrefix(kspInner, schurprefix);
899: /* Set DM for new solver */
900: KSPGetDM(jac->head->ksp, &dmInner);
901: KSPSetDM(kspInner, dmInner);
902: KSPSetDMActive(kspInner, PETSC_FALSE);
904: /* Defaults to PCKSP as preconditioner */
905: KSPGetPC(kspInner, &pcInner);
906: PCSetType(pcInner, PCKSP);
907: PCKSPSetKSP(pcInner, jac->head->ksp);
908: } else {
909: /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or
910: * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact,
911: * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for
912: * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make
913: * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used
914: * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */
915: KSPSetType(jac->head->ksp,KSPGMRES);
916: MatSchurComplementSetKSP(jac->schur,jac->head->ksp);
917: }
918: KSPSetOperators(jac->head->ksp,jac->mat[0],jac->pmat[0]);
919: KSPSetFromOptions(jac->head->ksp);
920: MatSetFromOptions(jac->schur);
922: PetscObjectTypeCompare((PetscObject)jac->schur, MATSCHURCOMPLEMENT, &flg);
923: if (flg) { /* Need to do this otherwise PCSetUp_KSP will overwrite the amat of jac->head->ksp */
924: KSP kspInner;
925: PC pcInner;
927: MatSchurComplementGetKSP(jac->schur, &kspInner);
928: KSPGetPC(kspInner, &pcInner);
929: PetscObjectTypeCompare((PetscObject)pcInner, PCKSP, &flg);
930: if (flg) {
931: KSP ksp;
933: PCKSPGetKSP(pcInner, &ksp);
934: if (ksp == jac->head->ksp) {
935: PCSetUseAmat(pcInner, PETSC_TRUE);
936: }
937: }
938: }
939: PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);
940: PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
941: if (flg) {
942: DM dmInner;
944: PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
945: KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper);
946: KSPSetErrorIfNotConverged(jac->kspupper,pc->erroriffailure);
947: KSPSetOptionsPrefix(jac->kspupper, schurprefix);
948: PetscObjectIncrementTabLevel((PetscObject)jac->kspupper, (PetscObject) pc, 1);
949: PetscObjectIncrementTabLevel((PetscObject)jac->kspupper->pc, (PetscObject) pc, 1);
950: KSPGetDM(jac->head->ksp, &dmInner);
951: KSPSetDM(jac->kspupper, dmInner);
952: KSPSetDMActive(jac->kspupper, PETSC_FALSE);
953: KSPSetFromOptions(jac->kspupper);
954: KSPSetOperators(jac->kspupper,jac->mat[0],jac->pmat[0]);
955: VecDuplicate(jac->head->x, &jac->head->z);
956: } else {
957: jac->kspupper = jac->head->ksp;
958: PetscObjectReference((PetscObject) jac->head->ksp);
959: }
961: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
962: MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);
963: }
964: KSPCreate(PetscObjectComm((PetscObject)pc),&jac->kspschur);
965: KSPSetErrorIfNotConverged(jac->kspschur,pc->erroriffailure);
966: PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);
967: PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);
968: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) {
969: PC pcschur;
970: KSPGetPC(jac->kspschur,&pcschur);
971: PCSetType(pcschur,PCNONE);
972: /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */
973: } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) {
974: MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user);
975: }
976: KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));
977: KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix);
978: KSPSetOptionsPrefix(jac->kspschur, Dprefix);
979: /* propagate DM */
980: {
981: DM sdm;
982: KSPGetDM(jac->head->next->ksp, &sdm);
983: if (sdm) {
984: KSPSetDM(jac->kspschur, sdm);
985: KSPSetDMActive(jac->kspschur, PETSC_FALSE);
986: }
987: }
988: /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
989: /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */
990: KSPSetFromOptions(jac->kspschur);
991: }
992: MatAssemblyBegin(jac->schur,MAT_FINAL_ASSEMBLY);
993: MatAssemblyEnd(jac->schur,MAT_FINAL_ASSEMBLY);
995: /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */
996: PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_L",ilink->splitname);
997: PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);
998: if (!LSC_L) {PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);}
999: if (LSC_L) {PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L);}
1000: PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_Lp",ilink->splitname);
1001: PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);
1002: if (!LSC_L) {PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);}
1003: if (LSC_L) {PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L);}
1004: } else if (jac->type == PC_COMPOSITE_GKB) {
1005: IS ccis;
1006: PetscInt rstart,rend;
1008: if (nsplit != 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"To use GKB preconditioner you must have exactly 2 fields");
1010: ilink = jac->head;
1012: /* When extracting off-diagonal submatrices, we take complements from this range */
1013: MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);
1015: ISComplement(ilink->is_col,rstart,rend,&ccis);
1016: if (jac->offdiag_use_amat) {
1017: MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
1018: } else {
1019: MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
1020: }
1021: ISDestroy(&ccis);
1022: /* Create work vectors for GKB algorithm */
1023: VecDuplicate(ilink->x,&jac->u);
1024: VecDuplicate(ilink->x,&jac->Hu);
1025: VecDuplicate(ilink->x,&jac->w2);
1026: ilink = ilink->next;
1027: ISComplement(ilink->is_col,rstart,rend,&ccis);
1028: if (jac->offdiag_use_amat) {
1029: MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
1030: } else {
1031: MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
1032: }
1033: ISDestroy(&ccis);
1034: /* Create work vectors for GKB algorithm */
1035: VecDuplicate(ilink->x,&jac->v);
1036: VecDuplicate(ilink->x,&jac->d);
1037: VecDuplicate(ilink->x,&jac->w1);
1038: MatGolubKahanComputeExplicitOperator(jac->mat[0],jac->B,jac->C,&jac->H,jac->gkbnu);
1039: PetscCalloc1(jac->gkbdelay,&jac->vecz);
1041: ilink = jac->head;
1042: KSPSetOperators(ilink->ksp,jac->H,jac->H);
1043: if (!jac->suboptionsset) {KSPSetFromOptions(ilink->ksp);}
1044: /* Create gkb_monitor context */
1045: if (jac->gkbmonitor) {
1046: PetscInt tablevel;
1047: PetscViewerCreate(PETSC_COMM_WORLD,&jac->gkbviewer);
1048: PetscViewerSetType(jac->gkbviewer,PETSCVIEWERASCII);
1049: PetscObjectGetTabLevel((PetscObject)ilink->ksp,&tablevel);
1050: PetscViewerASCIISetTab(jac->gkbviewer,tablevel);
1051: PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)ilink->ksp,1);
1052: }
1053: } else {
1054: /* set up the individual splits' PCs */
1055: i = 0;
1056: ilink = jac->head;
1057: while (ilink) {
1058: KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i]);
1059: /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
1060: if (!jac->suboptionsset) {KSPSetFromOptions(ilink->ksp);}
1061: i++;
1062: ilink = ilink->next;
1063: }
1064: }
1066: jac->suboptionsset = PETSC_TRUE;
1067: return(0);
1068: }
1070: #define FieldSplitSplitSolveAdd(ilink,xx,yy) \
1071: (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
1072: VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
1073: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\
1074: KSPSolve(ilink->ksp,ilink->x,ilink->y) || \
1075: KSPCheckSolve(ilink->ksp,pc,ilink->y) || \
1076: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\
1077: VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \
1078: VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE))
1080: static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y)
1081: {
1082: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1083: PetscErrorCode ierr;
1084: PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
1085: KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;
1088: switch (jac->schurfactorization) {
1089: case PC_FIELDSPLIT_SCHUR_FACT_DIAG:
1090: /* [A00 0; 0 -S], positive definite, suitable for MINRES */
1091: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1092: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
1093: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1094: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1095: KSPSolve(kspA,ilinkA->x,ilinkA->y);
1096: KSPCheckSolve(kspA,pc,ilinkA->y);
1097: PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1098: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1099: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
1100: PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1101: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
1102: KSPCheckSolve(jac->kspschur,pc,ilinkD->y);
1103: PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1104: VecScale(ilinkD->y,jac->schurscale);
1105: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1106: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1107: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1108: break;
1109: case PC_FIELDSPLIT_SCHUR_FACT_LOWER:
1110: /* [A00 0; A10 S], suitable for left preconditioning */
1111: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1112: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1113: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1114: KSPSolve(kspA,ilinkA->x,ilinkA->y);
1115: KSPCheckSolve(kspA,pc,ilinkA->y);
1116: PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1117: MatMult(jac->C,ilinkA->y,ilinkD->x);
1118: VecScale(ilinkD->x,-1.);
1119: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
1120: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1121: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
1122: PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1123: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
1124: KSPCheckSolve(jac->kspschur,pc,ilinkD->y);
1125: PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1126: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1127: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1128: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1129: break;
1130: case PC_FIELDSPLIT_SCHUR_FACT_UPPER:
1131: /* [A00 A01; 0 S], suitable for right preconditioning */
1132: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
1133: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
1134: PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1135: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
1136: KSPCheckSolve(jac->kspschur,pc,ilinkD->y);
1137: PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL); MatMult(jac->B,ilinkD->y,ilinkA->x);
1138: VecScale(ilinkA->x,-1.);
1139: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);
1140: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1141: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);
1142: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1143: KSPSolve(kspA,ilinkA->x,ilinkA->y);
1144: KSPCheckSolve(kspA,pc,ilinkA->y);
1145: PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1146: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1147: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1148: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1149: break;
1150: case PC_FIELDSPLIT_SCHUR_FACT_FULL:
1151: /* [1 0; A10 A00^{-1} 1] [A00 0; 0 S] [1 A00^{-1}A01; 0 1] */
1152: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1153: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1154: PetscLogEventBegin(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);
1155: KSPSolve(kspLower,ilinkA->x,ilinkA->y);
1156: KSPCheckSolve(kspLower,pc,ilinkA->y);
1157: PetscLogEventEnd(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);
1158: MatMult(jac->C,ilinkA->y,ilinkD->x);
1159: VecScale(ilinkD->x,-1.0);
1160: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
1161: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
1163: PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1164: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
1165: KSPCheckSolve(jac->kspschur,pc,ilinkD->y);
1166: PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1167: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1169: if (kspUpper == kspA) {
1170: MatMult(jac->B,ilinkD->y,ilinkA->y);
1171: VecAXPY(ilinkA->x,-1.0,ilinkA->y);
1172: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1173: KSPSolve(kspA,ilinkA->x,ilinkA->y);
1174: KSPCheckSolve(kspA,pc,ilinkA->y);
1175: PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1176: } else {
1177: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1178: KSPSolve(kspA,ilinkA->x,ilinkA->y);
1179: KSPCheckSolve(kspA,pc,ilinkA->y);
1180: MatMult(jac->B,ilinkD->y,ilinkA->x);
1181: PetscLogEventBegin(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);
1182: KSPSolve(kspUpper,ilinkA->x,ilinkA->z);
1183: KSPCheckSolve(kspUpper,pc,ilinkA->z);
1184: PetscLogEventEnd(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);
1185: VecAXPY(ilinkA->y,-1.0,ilinkA->z);
1186: }
1187: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1188: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1189: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1190: }
1191: return(0);
1192: }
1194: static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y)
1195: {
1196: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1197: PetscErrorCode ierr;
1198: PC_FieldSplitLink ilink = jac->head;
1199: PetscInt cnt,bs;
1202: if (jac->type == PC_COMPOSITE_ADDITIVE) {
1203: if (jac->defaultsplit) {
1204: VecGetBlockSize(x,&bs);
1205: 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);
1206: VecGetBlockSize(y,&bs);
1207: 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);
1208: VecStrideGatherAll(x,jac->x,INSERT_VALUES);
1209: while (ilink) {
1210: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1211: KSPSolve(ilink->ksp,ilink->x,ilink->y);
1212: KSPCheckSolve(ilink->ksp,pc,ilink->y);
1213: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1214: ilink = ilink->next;
1215: }
1216: VecStrideScatterAll(jac->y,y,INSERT_VALUES);
1217: } else {
1218: VecSet(y,0.0);
1219: while (ilink) {
1220: FieldSplitSplitSolveAdd(ilink,x,y);
1221: ilink = ilink->next;
1222: }
1223: }
1224: } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) {
1225: VecSet(y,0.0);
1226: /* solve on first block for first block variables */
1227: VecScatterBegin(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);
1228: VecScatterEnd(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);
1229: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1230: KSPSolve(ilink->ksp,ilink->x,ilink->y);
1231: KSPCheckSolve(ilink->ksp,pc,ilink->y);
1232: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1233: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1234: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1236: /* compute the residual only onto second block variables using first block variables */
1237: MatMult(jac->Afield[1],ilink->y,ilink->next->x);
1238: ilink = ilink->next;
1239: VecScale(ilink->x,-1.0);
1240: VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1241: VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1243: /* solve on second block variables */
1244: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1245: KSPSolve(ilink->ksp,ilink->x,ilink->y);
1246: KSPCheckSolve(ilink->ksp,pc,ilink->y);
1247: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1248: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1249: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1250: } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1251: if (!jac->w1) {
1252: VecDuplicate(x,&jac->w1);
1253: VecDuplicate(x,&jac->w2);
1254: }
1255: VecSet(y,0.0);
1256: FieldSplitSplitSolveAdd(ilink,x,y);
1257: cnt = 1;
1258: while (ilink->next) {
1259: ilink = ilink->next;
1260: /* compute the residual only over the part of the vector needed */
1261: MatMult(jac->Afield[cnt++],y,ilink->x);
1262: VecScale(ilink->x,-1.0);
1263: VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1264: VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1265: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1266: KSPSolve(ilink->ksp,ilink->x,ilink->y);
1267: KSPCheckSolve(ilink->ksp,pc,ilink->y);
1268: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1269: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1270: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1271: }
1272: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1273: cnt -= 2;
1274: while (ilink->previous) {
1275: ilink = ilink->previous;
1276: /* compute the residual only over the part of the vector needed */
1277: MatMult(jac->Afield[cnt--],y,ilink->x);
1278: VecScale(ilink->x,-1.0);
1279: VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1280: VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1281: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1282: KSPSolve(ilink->ksp,ilink->x,ilink->y);
1283: KSPCheckSolve(ilink->ksp,pc,ilink->y);
1284: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1285: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1286: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1287: }
1288: }
1289: } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type);
1290: return(0);
1291: }
1294: static PetscErrorCode PCApply_FieldSplit_GKB(PC pc,Vec x,Vec y)
1295: {
1296: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1297: PetscErrorCode ierr;
1298: PC_FieldSplitLink ilinkA = jac->head,ilinkD = ilinkA->next;
1299: KSP ksp = ilinkA->ksp;
1300: Vec u,v,Hu,d,work1,work2;
1301: PetscScalar alpha,z,nrmz2,*vecz;
1302: PetscReal lowbnd,nu,beta;
1303: PetscInt j,iterGKB;
1306: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1307: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
1308: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1309: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
1311: u = jac->u;
1312: v = jac->v;
1313: Hu = jac->Hu;
1314: d = jac->d;
1315: work1 = jac->w1;
1316: work2 = jac->w2;
1317: vecz = jac->vecz;
1319: /* Change RHS to comply with matrix regularization H = A + nu*B*B' */
1320: /* Add q = q + nu*B*b */
1321: if (jac->gkbnu) {
1322: nu = jac->gkbnu;
1323: VecScale(ilinkD->x,jac->gkbnu);
1324: MatMultAdd(jac->B,ilinkD->x,ilinkA->x,ilinkA->x); /* q = q + nu*B*b */
1325: } else {
1326: /* Situation when no augmented Lagrangian is used. Then we set inner */
1327: /* matrix N = I in [Ar13], and thus nu = 1. */
1328: nu = 1;
1329: }
1331: /* Transform rhs from [q,tilde{b}] to [0,b] */
1332: PetscLogEventBegin(ilinkA->event,ksp,ilinkA->x,ilinkA->y,NULL);
1333: KSPSolve(ksp,ilinkA->x,ilinkA->y);
1334: KSPCheckSolve(ksp,pc,ilinkA->y);
1335: PetscLogEventEnd(ilinkA->event,ksp,ilinkA->x,ilinkA->y,NULL);
1336: MatMultHermitianTranspose(jac->B,ilinkA->y,work1);
1337: VecAXPBY(work1,1.0/nu,-1.0,ilinkD->x); /* c = b - B'*x */
1339: /* First step of algorithm */
1340: VecNorm(work1,NORM_2,&beta); /* beta = sqrt(nu*c'*c)*/
1341: KSPCheckDot(ksp,beta);
1342: beta = PetscSqrtScalar(nu)*beta;
1343: VecAXPBY(v,nu/beta,0.0,work1); /* v = nu/beta *c */
1344: MatMult(jac->B,v,work2); /* u = H^{-1}*B*v */
1345: PetscLogEventBegin(ilinkA->event,ksp,work2,u,NULL);
1346: KSPSolve(ksp,work2,u);
1347: KSPCheckSolve(ksp,pc,u);
1348: PetscLogEventEnd(ilinkA->event,ksp,work2,u,NULL);
1349: MatMult(jac->H,u,Hu); /* alpha = u'*H*u */
1350: VecDot(Hu,u,&alpha);
1351: KSPCheckDot(ksp,alpha);
1352: if (PetscRealPart(alpha) <= 0.0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_NOT_CONVERGED,"GKB preconditioner diverged, H is not positive definite");
1353: alpha = PetscSqrtScalar(PetscAbsScalar(alpha));
1354: VecScale(u,1.0/alpha);
1355: VecAXPBY(d,1.0/alpha,0.0,v); /* v = nu/beta *c */
1357: z = beta/alpha;
1358: vecz[1] = z;
1360: /* Computation of first iterate x(1) and p(1) */
1361: VecAXPY(ilinkA->y,z,u);
1362: VecCopy(d,ilinkD->y);
1363: VecScale(ilinkD->y,-z);
1365: iterGKB = 1; lowbnd = 2*jac->gkbtol;
1366: if (jac->gkbmonitor) {
1367: PetscViewerASCIIPrintf(jac->gkbviewer,"%3D GKB Lower bound estimate %14.12e\n",iterGKB,lowbnd);
1368: }
1370: while (iterGKB < jac->gkbmaxit && lowbnd > jac->gkbtol) {
1371: iterGKB += 1;
1372: MatMultHermitianTranspose(jac->B,u,work1); /* v <- nu*(B'*u-alpha/nu*v) */
1373: VecAXPBY(v,nu,-alpha,work1);
1374: VecNorm(v,NORM_2,&beta); /* beta = sqrt(nu)*v'*v */
1375: beta = beta/PetscSqrtScalar(nu);
1376: VecScale(v,1.0/beta);
1377: MatMult(jac->B,v,work2); /* u <- H^{-1}*(B*v-beta*H*u) */
1378: MatMult(jac->H,u,Hu);
1379: VecAXPY(work2,-beta,Hu);
1380: PetscLogEventBegin(ilinkA->event,ksp,work2,u,NULL);
1381: KSPSolve(ksp,work2,u);
1382: KSPCheckSolve(ksp,pc,u);
1383: PetscLogEventEnd(ilinkA->event,ksp,work2,u,NULL);
1384: MatMult(jac->H,u,Hu); /* alpha = u'*H*u */
1385: VecDot(Hu,u,&alpha);
1386: KSPCheckDot(ksp,alpha);
1387: if (PetscRealPart(alpha) <= 0.0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_NOT_CONVERGED,"GKB preconditioner diverged, H is not positive definite");
1388: alpha = PetscSqrtScalar(PetscAbsScalar(alpha));
1389: VecScale(u,1.0/alpha);
1391: z = -beta/alpha*z; /* z <- beta/alpha*z */
1392: vecz[0] = z;
1394: /* Computation of new iterate x(i+1) and p(i+1) */
1395: VecAXPBY(d,1.0/alpha,-beta/alpha,v); /* d = (v-beta*d)/alpha */
1396: VecAXPY(ilinkA->y,z,u); /* r = r + z*u */
1397: VecAXPY(ilinkD->y,-z,d); /* p = p - z*d */
1398: MatMult(jac->H,ilinkA->y,Hu); /* ||u||_H = u'*H*u */
1399: VecDot(Hu,ilinkA->y,&nrmz2);
1401: /* Compute Lower Bound estimate */
1402: if (iterGKB > jac->gkbdelay) {
1403: lowbnd = 0.0;
1404: for (j=0; j<jac->gkbdelay; j++) {
1405: lowbnd += PetscAbsScalar(vecz[j]*vecz[j]);
1406: }
1407: lowbnd = PetscSqrtScalar(lowbnd/PetscAbsScalar(nrmz2));
1408: }
1410: for (j=0; j<jac->gkbdelay-1; j++) {
1411: vecz[jac->gkbdelay-j-1] = vecz[jac->gkbdelay-j-2];
1412: }
1413: if (jac->gkbmonitor) {
1414: PetscViewerASCIIPrintf(jac->gkbviewer,"%3D GKB Lower bound estimate %14.12e\n",iterGKB,lowbnd);
1415: }
1416: }
1418: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1419: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1420: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1421: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1423: return(0);
1424: }
1427: #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \
1428: (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
1429: VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
1430: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \
1431: KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \
1432: KSPCheckSolve(ilink->ksp,pc,ilink->x) || \
1433: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \
1434: VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \
1435: VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE))
1437: static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y)
1438: {
1439: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1440: PetscErrorCode ierr;
1441: PC_FieldSplitLink ilink = jac->head;
1442: PetscInt bs;
1445: if (jac->type == PC_COMPOSITE_ADDITIVE) {
1446: if (jac->defaultsplit) {
1447: VecGetBlockSize(x,&bs);
1448: 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);
1449: VecGetBlockSize(y,&bs);
1450: 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);
1451: VecStrideGatherAll(x,jac->x,INSERT_VALUES);
1452: while (ilink) {
1453: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1454: KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);
1455: KSPCheckSolve(ilink->ksp,pc,ilink->y);
1456: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1457: ilink = ilink->next;
1458: }
1459: VecStrideScatterAll(jac->y,y,INSERT_VALUES);
1460: } else {
1461: VecSet(y,0.0);
1462: while (ilink) {
1463: FieldSplitSplitSolveAddTranspose(ilink,x,y);
1464: ilink = ilink->next;
1465: }
1466: }
1467: } else {
1468: if (!jac->w1) {
1469: VecDuplicate(x,&jac->w1);
1470: VecDuplicate(x,&jac->w2);
1471: }
1472: VecSet(y,0.0);
1473: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1474: FieldSplitSplitSolveAddTranspose(ilink,x,y);
1475: while (ilink->next) {
1476: ilink = ilink->next;
1477: MatMultTranspose(pc->mat,y,jac->w1);
1478: VecWAXPY(jac->w2,-1.0,jac->w1,x);
1479: FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1480: }
1481: while (ilink->previous) {
1482: ilink = ilink->previous;
1483: MatMultTranspose(pc->mat,y,jac->w1);
1484: VecWAXPY(jac->w2,-1.0,jac->w1,x);
1485: FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1486: }
1487: } else {
1488: while (ilink->next) { /* get to last entry in linked list */
1489: ilink = ilink->next;
1490: }
1491: FieldSplitSplitSolveAddTranspose(ilink,x,y);
1492: while (ilink->previous) {
1493: ilink = ilink->previous;
1494: MatMultTranspose(pc->mat,y,jac->w1);
1495: VecWAXPY(jac->w2,-1.0,jac->w1,x);
1496: FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1497: }
1498: }
1499: }
1500: return(0);
1501: }
1503: static PetscErrorCode PCReset_FieldSplit(PC pc)
1504: {
1505: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1506: PetscErrorCode ierr;
1507: PC_FieldSplitLink ilink = jac->head,next;
1510: while (ilink) {
1511: KSPDestroy(&ilink->ksp);
1512: VecDestroy(&ilink->x);
1513: VecDestroy(&ilink->y);
1514: VecDestroy(&ilink->z);
1515: VecScatterDestroy(&ilink->sctx);
1516: ISDestroy(&ilink->is);
1517: ISDestroy(&ilink->is_col);
1518: PetscFree(ilink->splitname);
1519: PetscFree(ilink->fields);
1520: PetscFree(ilink->fields_col);
1521: next = ilink->next;
1522: PetscFree(ilink);
1523: ilink = next;
1524: }
1525: jac->head = NULL;
1526: PetscFree2(jac->x,jac->y);
1527: if (jac->mat && jac->mat != jac->pmat) {
1528: MatDestroyMatrices(jac->nsplits,&jac->mat);
1529: } else if (jac->mat) {
1530: jac->mat = NULL;
1531: }
1532: if (jac->pmat) {MatDestroyMatrices(jac->nsplits,&jac->pmat);}
1533: if (jac->Afield) {MatDestroyMatrices(jac->nsplits,&jac->Afield);}
1534: jac->nsplits = 0;
1535: VecDestroy(&jac->w1);
1536: VecDestroy(&jac->w2);
1537: MatDestroy(&jac->schur);
1538: MatDestroy(&jac->schurp);
1539: MatDestroy(&jac->schur_user);
1540: KSPDestroy(&jac->kspschur);
1541: KSPDestroy(&jac->kspupper);
1542: MatDestroy(&jac->B);
1543: MatDestroy(&jac->C);
1544: MatDestroy(&jac->H);
1545: VecDestroy(&jac->u);
1546: VecDestroy(&jac->v);
1547: VecDestroy(&jac->Hu);
1548: VecDestroy(&jac->d);
1549: PetscFree(jac->vecz);
1550: PetscViewerDestroy(&jac->gkbviewer);
1551: jac->isrestrict = PETSC_FALSE;
1552: return(0);
1553: }
1555: static PetscErrorCode PCDestroy_FieldSplit(PC pc)
1556: {
1557: PetscErrorCode ierr;
1560: PCReset_FieldSplit(pc);
1561: PetscFree(pc->data);
1562: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurGetSubKSP_C",NULL);
1563: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL);
1564: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL);
1565: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL);
1566: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL);
1567: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL);
1568: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL);
1569: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL);
1570: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL);
1571: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",NULL);
1572: return(0);
1573: }
1575: static PetscErrorCode PCSetFromOptions_FieldSplit(PetscOptionItems *PetscOptionsObject,PC pc)
1576: {
1577: PetscErrorCode ierr;
1578: PetscInt bs;
1579: PetscBool flg;
1580: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1581: PCCompositeType ctype;
1584: PetscOptionsHead(PetscOptionsObject,"FieldSplit options");
1585: PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL);
1586: PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);
1587: if (flg) {
1588: PCFieldSplitSetBlockSize(pc,bs);
1589: }
1590: jac->diag_use_amat = pc->useAmat;
1591: 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);
1592: jac->offdiag_use_amat = pc->useAmat;
1593: 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);
1594: PetscOptionsBool("-pc_fieldsplit_detect_saddle_point","Form 2-way split by detecting zero diagonal entries", "PCFieldSplitSetDetectSaddlePoint",jac->detect,&jac->detect,NULL);
1595: PCFieldSplitSetDetectSaddlePoint(pc,jac->detect); /* Sets split type and Schur PC type */
1596: PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);
1597: if (flg) {
1598: PCFieldSplitSetType(pc,ctype);
1599: }
1600: /* Only setup fields once */
1601: if ((jac->bs > 0) && (jac->nsplits == 0)) {
1602: /* only allow user to set fields from command line if bs is already known.
1603: otherwise user can set them in PCFieldSplitSetDefaults() */
1604: PCFieldSplitSetRuntimeSplits_Private(pc);
1605: if (jac->splitdefined) {PetscInfo(pc,"Splits defined using the options database\n");}
1606: }
1607: if (jac->type == PC_COMPOSITE_SCHUR) {
1608: PetscOptionsGetEnum(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);
1609: if (flg) {PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");}
1610: 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);
1611: PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSetSchurPre",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL);
1612: PetscOptionsScalar("-pc_fieldsplit_schur_scale","Scale Schur complement","PCFieldSplitSetSchurScale",jac->schurscale,&jac->schurscale,NULL);
1613: } else if (jac->type == PC_COMPOSITE_GKB) {
1614: PetscOptionsReal("-pc_fieldsplit_gkb_tol","The tolerance for the lower bound stopping criterion","PCFieldSplitGKBTol",jac->gkbtol,&jac->gkbtol,NULL);
1615: PetscOptionsInt("-pc_fieldsplit_gkb_delay","The delay value for lower bound criterion","PCFieldSplitGKBDelay",jac->gkbdelay,&jac->gkbdelay,NULL);
1616: PetscOptionsReal("-pc_fieldsplit_gkb_nu","Parameter in augmented Lagrangian approach","PCFieldSplitGKBNu",jac->gkbnu,&jac->gkbnu,NULL);
1617: if (jac->gkbnu < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nu cannot be less than 0: value %f",jac->gkbnu);
1618: PetscOptionsInt("-pc_fieldsplit_gkb_maxit","Maximum allowed number of iterations","PCFieldSplitGKBMaxit",jac->gkbmaxit,&jac->gkbmaxit,NULL);
1619: PetscOptionsBool("-pc_fieldsplit_gkb_monitor","Prints number of GKB iterations and error","PCFieldSplitGKB",jac->gkbmonitor,&jac->gkbmonitor,NULL);
1620: }
1621: PetscOptionsTail();
1622: return(0);
1623: }
1625: /*------------------------------------------------------------------------------------*/
1627: static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1628: {
1629: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1630: PetscErrorCode ierr;
1631: PC_FieldSplitLink ilink,next = jac->head;
1632: char prefix[128];
1633: PetscInt i;
1636: if (jac->splitdefined) {
1637: PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1638: return(0);
1639: }
1640: for (i=0; i<n; i++) {
1641: if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs);
1642: if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]);
1643: }
1644: PetscNew(&ilink);
1645: if (splitname) {
1646: PetscStrallocpy(splitname,&ilink->splitname);
1647: } else {
1648: PetscMalloc1(3,&ilink->splitname);
1649: PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);
1650: }
1651: ilink->event = jac->nsplits < 5 ? KSP_Solve_FS_0 + jac->nsplits : KSP_Solve_FS_0 + 4; /* Any split great than 4 gets logged in the 4th split */
1652: PetscMalloc1(n,&ilink->fields);
1653: PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));
1654: PetscMalloc1(n,&ilink->fields_col);
1655: PetscMemcpy(ilink->fields_col,fields_col,n*sizeof(PetscInt));
1657: ilink->nfields = n;
1658: ilink->next = NULL;
1659: KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1660: KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);
1661: PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1662: KSPSetType(ilink->ksp,KSPPREONLY);
1663: PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);
1665: PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);
1666: KSPSetOptionsPrefix(ilink->ksp,prefix);
1668: if (!next) {
1669: jac->head = ilink;
1670: ilink->previous = NULL;
1671: } else {
1672: while (next->next) {
1673: next = next->next;
1674: }
1675: next->next = ilink;
1676: ilink->previous = next;
1677: }
1678: jac->nsplits++;
1679: return(0);
1680: }
1682: static PetscErrorCode PCFieldSplitSchurGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp)
1683: {
1684: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1688: *subksp = NULL;
1689: if (n) *n = 0;
1690: if (jac->type == PC_COMPOSITE_SCHUR) {
1691: PetscInt nn;
1693: if (!jac->schur) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitSchurGetSubKSP()");
1694: if (jac->nsplits != 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Unexpected number of splits %D != 2",jac->nsplits);
1695: nn = jac->nsplits + (jac->kspupper != jac->head->ksp ? 1 : 0);
1696: PetscMalloc1(nn,subksp);
1697: (*subksp)[0] = jac->head->ksp;
1698: (*subksp)[1] = jac->kspschur;
1699: if (jac->kspupper != jac->head->ksp) (*subksp)[2] = jac->kspupper;
1700: if (n) *n = nn;
1701: }
1702: return(0);
1703: }
1705: static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp)
1706: {
1707: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1711: if (!jac->schur) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitGetSubKSP()");
1712: PetscMalloc1(jac->nsplits,subksp);
1713: MatSchurComplementGetKSP(jac->schur,*subksp);
1715: (*subksp)[1] = jac->kspschur;
1716: if (n) *n = jac->nsplits;
1717: return(0);
1718: }
1720: static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp)
1721: {
1722: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1723: PetscErrorCode ierr;
1724: PetscInt cnt = 0;
1725: PC_FieldSplitLink ilink = jac->head;
1728: PetscMalloc1(jac->nsplits,subksp);
1729: while (ilink) {
1730: (*subksp)[cnt++] = ilink->ksp;
1731: ilink = ilink->next;
1732: }
1733: 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);
1734: if (n) *n = jac->nsplits;
1735: return(0);
1736: }
1738: /*@C
1739: PCFieldSplitRestrictIS - Restricts the fieldsplit ISs to be within a given IS.
1741: Input Parameters:
1742: + pc - the preconditioner context
1743: + is - the index set that defines the indices to which the fieldsplit is to be restricted
1745: Level: advanced
1747: @*/
1748: PetscErrorCode PCFieldSplitRestrictIS(PC pc,IS isy)
1749: {
1755: PetscTryMethod(pc,"PCFieldSplitRestrictIS_C",(PC,IS),(pc,isy));
1756: return(0);
1757: }
1760: static PetscErrorCode PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy)
1761: {
1762: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1763: PetscErrorCode ierr;
1764: PC_FieldSplitLink ilink = jac->head, next;
1765: PetscInt localsize,size,sizez,i;
1766: const PetscInt *ind, *indz;
1767: PetscInt *indc, *indcz;
1768: PetscBool flg;
1771: ISGetLocalSize(isy,&localsize);
1772: MPI_Scan(&localsize,&size,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)isy));
1773: size -= localsize;
1774: while(ilink) {
1775: IS isrl,isr;
1776: PC subpc;
1777: ISEmbed(ilink->is, isy, PETSC_TRUE, &isrl);
1778: ISGetLocalSize(isrl,&localsize);
1779: PetscMalloc1(localsize,&indc);
1780: ISGetIndices(isrl,&ind);
1781: PetscMemcpy(indc,ind,localsize*sizeof(PetscInt));
1782: ISRestoreIndices(isrl,&ind);
1783: ISDestroy(&isrl);
1784: for (i=0; i<localsize; i++) *(indc+i) += size;
1785: ISCreateGeneral(PetscObjectComm((PetscObject)isy),localsize,indc,PETSC_OWN_POINTER,&isr);
1786: PetscObjectReference((PetscObject)isr);
1787: ISDestroy(&ilink->is);
1788: ilink->is = isr;
1789: PetscObjectReference((PetscObject)isr);
1790: ISDestroy(&ilink->is_col);
1791: ilink->is_col = isr;
1792: ISDestroy(&isr);
1793: KSPGetPC(ilink->ksp, &subpc);
1794: PetscObjectTypeCompare((PetscObject)subpc,PCFIELDSPLIT,&flg);
1795: if(flg) {
1796: IS iszl,isz;
1797: MPI_Comm comm;
1798: ISGetLocalSize(ilink->is,&localsize);
1799: comm = PetscObjectComm((PetscObject)ilink->is);
1800: ISEmbed(isy, ilink->is, PETSC_TRUE, &iszl);
1801: MPI_Scan(&localsize,&sizez,1,MPIU_INT,MPI_SUM,comm);
1802: sizez -= localsize;
1803: ISGetLocalSize(iszl,&localsize);
1804: PetscMalloc1(localsize,&indcz);
1805: ISGetIndices(iszl,&indz);
1806: PetscMemcpy(indcz,indz,localsize*sizeof(PetscInt));
1807: ISRestoreIndices(iszl,&indz);
1808: ISDestroy(&iszl);
1809: for (i=0; i<localsize; i++) *(indcz+i) += sizez;
1810: ISCreateGeneral(comm,localsize,indcz,PETSC_OWN_POINTER,&isz);
1811: PCFieldSplitRestrictIS(subpc,isz);
1812: ISDestroy(&isz);
1813: }
1814: next = ilink->next;
1815: ilink = next;
1816: }
1817: jac->isrestrict = PETSC_TRUE;
1818: return(0);
1819: }
1821: static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is)
1822: {
1823: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1824: PetscErrorCode ierr;
1825: PC_FieldSplitLink ilink, next = jac->head;
1826: char prefix[128];
1829: if (jac->splitdefined) {
1830: PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1831: return(0);
1832: }
1833: PetscNew(&ilink);
1834: if (splitname) {
1835: PetscStrallocpy(splitname,&ilink->splitname);
1836: } else {
1837: PetscMalloc1(8,&ilink->splitname);
1838: PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);
1839: }
1840: ilink->event = jac->nsplits < 5 ? KSP_Solve_FS_0 + jac->nsplits : KSP_Solve_FS_0 + 4; /* Any split great than 4 gets logged in the 4th split */
1841: PetscObjectReference((PetscObject)is);
1842: ISDestroy(&ilink->is);
1843: ilink->is = is;
1844: PetscObjectReference((PetscObject)is);
1845: ISDestroy(&ilink->is_col);
1846: ilink->is_col = is;
1847: ilink->next = NULL;
1848: KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1849: KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);
1850: PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1851: KSPSetType(ilink->ksp,KSPPREONLY);
1852: PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);
1854: PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);
1855: KSPSetOptionsPrefix(ilink->ksp,prefix);
1857: if (!next) {
1858: jac->head = ilink;
1859: ilink->previous = NULL;
1860: } else {
1861: while (next->next) {
1862: next = next->next;
1863: }
1864: next->next = ilink;
1865: ilink->previous = next;
1866: }
1867: jac->nsplits++;
1868: return(0);
1869: }
1871: /*@C
1872: PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner
1874: Logically Collective on PC
1876: Input Parameters:
1877: + pc - the preconditioner context
1878: . splitname - name of this split, if NULL the number of the split is used
1879: . n - the number of fields in this split
1880: - fields - the fields in this split
1882: Level: intermediate
1884: Notes:
1885: Use PCFieldSplitSetIS() to set a completely general set of indices as a field.
1887: The PCFieldSplitSetFields() is for defining fields as strided blocks. For example, if the block
1888: 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
1889: 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x....
1890: where the numbered entries indicate what is in the field.
1892: This function is called once per split (it creates a new split each time). Solve options
1893: for this split will be available under the prefix -fieldsplit_SPLITNAME_.
1895: Developer Note: This routine does not actually create the IS representing the split, that is delayed
1896: until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be
1897: available when this routine is called.
1899: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS()
1901: @*/
1902: PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1903: {
1909: if (n < 1) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname);
1911: PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt*,const PetscInt*),(pc,splitname,n,fields,fields_col));
1912: return(0);
1913: }
1915: /*@
1916: PCFieldSplitSetDiagUseAmat - set flag indicating whether to extract diagonal blocks from Amat (rather than Pmat)
1918: Logically Collective on PC
1920: Input Parameters:
1921: + pc - the preconditioner object
1922: - flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from
1924: Options Database:
1925: . -pc_fieldsplit_diag_use_amat
1927: Level: intermediate
1929: .seealso: PCFieldSplitGetDiagUseAmat(), PCFieldSplitSetOffDiagUseAmat(), PCFIELDSPLIT
1931: @*/
1932: PetscErrorCode PCFieldSplitSetDiagUseAmat(PC pc,PetscBool flg)
1933: {
1934: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1935: PetscBool isfs;
1940: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1941: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1942: jac->diag_use_amat = flg;
1943: return(0);
1944: }
1946: /*@
1947: PCFieldSplitGetDiagUseAmat - get the flag indicating whether to extract diagonal blocks from Amat (rather than Pmat)
1949: Logically Collective on PC
1951: Input Parameters:
1952: . pc - the preconditioner object
1954: Output Parameters:
1955: . flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from
1958: Level: intermediate
1960: .seealso: PCFieldSplitSetDiagUseAmat(), PCFieldSplitGetOffDiagUseAmat(), PCFIELDSPLIT
1962: @*/
1963: PetscErrorCode PCFieldSplitGetDiagUseAmat(PC pc,PetscBool *flg)
1964: {
1965: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1966: PetscBool isfs;
1972: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1973: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1974: *flg = jac->diag_use_amat;
1975: return(0);
1976: }
1978: /*@
1979: PCFieldSplitSetOffDiagUseAmat - set flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat)
1981: Logically Collective on PC
1983: Input Parameters:
1984: + pc - the preconditioner object
1985: - flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from
1987: Options Database:
1988: . -pc_fieldsplit_off_diag_use_amat
1990: Level: intermediate
1992: .seealso: PCFieldSplitGetOffDiagUseAmat(), PCFieldSplitSetDiagUseAmat(), PCFIELDSPLIT
1994: @*/
1995: PetscErrorCode PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg)
1996: {
1997: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1998: PetscBool isfs;
2003: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2004: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
2005: jac->offdiag_use_amat = flg;
2006: return(0);
2007: }
2009: /*@
2010: PCFieldSplitGetOffDiagUseAmat - get the flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat)
2012: Logically Collective on PC
2014: Input Parameters:
2015: . pc - the preconditioner object
2017: Output Parameters:
2018: . flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from
2021: Level: intermediate
2023: .seealso: PCFieldSplitSetOffDiagUseAmat(), PCFieldSplitGetDiagUseAmat(), PCFIELDSPLIT
2025: @*/
2026: PetscErrorCode PCFieldSplitGetOffDiagUseAmat(PC pc,PetscBool *flg)
2027: {
2028: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2029: PetscBool isfs;
2035: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2036: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
2037: *flg = jac->offdiag_use_amat;
2038: return(0);
2039: }
2043: /*@C
2044: PCFieldSplitSetIS - Sets the exact elements for field
2046: Logically Collective on PC
2048: Input Parameters:
2049: + pc - the preconditioner context
2050: . splitname - name of this split, if NULL the number of the split is used
2051: - is - the index set that defines the vector elements in this field
2054: Notes:
2055: Use PCFieldSplitSetFields(), for fields defined by strided types.
2057: This function is called once per split (it creates a new split each time). Solve options
2058: for this split will be available under the prefix -fieldsplit_SPLITNAME_.
2060: Level: intermediate
2062: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize()
2064: @*/
2065: PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is)
2066: {
2073: PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));
2074: return(0);
2075: }
2077: /*@C
2078: PCFieldSplitGetIS - Retrieves the elements for a field as an IS
2080: Logically Collective on PC
2082: Input Parameters:
2083: + pc - the preconditioner context
2084: - splitname - name of this split
2086: Output Parameter:
2087: - is - the index set that defines the vector elements in this field, or NULL if the field is not found
2089: Level: intermediate
2091: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS()
2093: @*/
2094: PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is)
2095: {
2102: {
2103: PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;
2104: PC_FieldSplitLink ilink = jac->head;
2105: PetscBool found;
2107: *is = NULL;
2108: while (ilink) {
2109: PetscStrcmp(ilink->splitname, splitname, &found);
2110: if (found) {
2111: *is = ilink->is;
2112: break;
2113: }
2114: ilink = ilink->next;
2115: }
2116: }
2117: return(0);
2118: }
2120: /*@
2121: PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the
2122: fieldsplit preconditioner. If not set the matrix block size is used.
2124: Logically Collective on PC
2126: Input Parameters:
2127: + pc - the preconditioner context
2128: - bs - the block size
2130: Level: intermediate
2132: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields()
2134: @*/
2135: PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs)
2136: {
2142: PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));
2143: return(0);
2144: }
2146: /*@C
2147: PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits
2149: Collective on KSP
2151: Input Parameter:
2152: . pc - the preconditioner context
2154: Output Parameters:
2155: + n - the number of splits
2156: - subksp - the array of KSP contexts
2158: Note:
2159: After PCFieldSplitGetSubKSP() the array of KSPs is to be freed by the user with PetscFree()
2160: (not the KSP just the array that contains them).
2162: You must call PCSetUp() before calling PCFieldSplitGetSubKSP().
2164: If the fieldsplit is of type PC_COMPOSITE_SCHUR, it returns the KSP object used inside the
2165: Schur complement and the KSP object used to iterate over the Schur complement.
2166: To access all the KSP objects used in PC_COMPOSITE_SCHUR, use PCFieldSplitSchurGetSubKSP().
2168: If the fieldsplit is of type PC_COMPOSITE_GKB, it returns the KSP object used to solve the
2169: inner linear system defined by the matrix H in each loop.
2171: Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs.
2172: You can call PCFieldSplitGetSubKSP(pc,n,PETSC_NULL_KSP,ierr) to determine how large the
2173: KSP array must be.
2176: Level: advanced
2178: .seealso: PCFIELDSPLIT
2179: @*/
2180: PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
2181: {
2187: PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));
2188: return(0);
2189: }
2191: /*@C
2192: PCFieldSplitSchurGetSubKSP - Gets the KSP contexts used inside the Schur complement based PCFIELDSPLIT
2194: Collective on KSP
2196: Input Parameter:
2197: . pc - the preconditioner context
2199: Output Parameters:
2200: + n - the number of splits
2201: - subksp - the array of KSP contexts
2203: Note:
2204: After PCFieldSplitSchurGetSubKSP() the array of KSPs is to be freed by the user with PetscFree()
2205: (not the KSP just the array that contains them).
2207: You must call PCSetUp() before calling PCFieldSplitSchurGetSubKSP().
2209: If the fieldsplit type is of type PC_COMPOSITE_SCHUR, it returns (in order)
2210: - the KSP used for the (1,1) block
2211: - the KSP used for the Schur complement (not the one used for the interior Schur solver)
2212: - the KSP used for the (1,1) block in the upper triangular factor (if different from that of the (1,1) block).
2214: It returns a null array if the fieldsplit is not of type PC_COMPOSITE_SCHUR; in this case, you should use PCFieldSplitGetSubKSP().
2216: Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs.
2217: You can call PCFieldSplitSchurGetSubKSP(pc,n,PETSC_NULL_KSP,ierr) to determine how large the
2218: KSP array must be.
2220: Level: advanced
2222: .seealso: PCFIELDSPLIT
2223: @*/
2224: PetscErrorCode PCFieldSplitSchurGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
2225: {
2231: PetscUseMethod(pc,"PCFieldSplitSchurGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));
2232: return(0);
2233: }
2235: /*@
2236: PCFieldSplitSetSchurPre - Indicates what operator is used to construct the preconditioner for the Schur complement.
2237: A11 matrix. Otherwise no preconditioner is used.
2239: Collective on PC
2241: Input Parameters:
2242: + pc - the preconditioner context
2243: . ptype - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11 (default), PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_SCHUR_PRE_USER
2244: PC_FIELDSPLIT_SCHUR_PRE_SELFP, and PC_FIELDSPLIT_SCHUR_PRE_FULL
2245: - userpre - matrix to use for preconditioning, or NULL
2247: Options Database:
2248: . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11. See notes for meaning of various arguments
2250: Notes:
2251: $ If ptype is
2252: $ a11 then the preconditioner for the Schur complement is generated from the block diagonal part of the preconditioner
2253: $ matrix associated with the Schur complement (i.e. A11), not the Schur complement matrix
2254: $ self the preconditioner for the Schur complement is generated from the symbolic representation of the Schur complement matrix:
2255: $ The only preconditioner that currently works with this symbolic respresentation matrix object is the PCLSC
2256: $ preconditioner
2257: $ user then the preconditioner for the Schur complement is generated from the user provided matrix (pre argument
2258: $ to this function).
2259: $ selfp then the preconditioning for the Schur complement is generated from an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01
2260: $ This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be
2261: $ lumped before extracting the diagonal using the additional option -fieldsplit_1_mat_schur_complement_ainv_type lump
2262: $ full then the preconditioner for the Schur complement is generated from the exact Schur complement matrix representation computed internally by PCFIELDSPLIT (this is expensive)
2263: $ useful mostly as a test that the Schur complement approach can work for your problem
2265: When solving a saddle point problem, where the A11 block is identically zero, using a11 as the ptype only makes sense
2266: with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and
2267: -fieldsplit_1_pc_type lsc which uses the least squares commutator to compute a preconditioner for the Schur complement.
2269: Level: intermediate
2271: .seealso: PCFieldSplitGetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType,
2272: MatSchurComplementSetAinvType(), PCLSC
2274: @*/
2275: PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
2276: {
2281: PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));
2282: return(0);
2283: }
2285: PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) {return PCFieldSplitSetSchurPre(pc,ptype,pre);} /* Deprecated name */
2287: /*@
2288: PCFieldSplitGetSchurPre - For Schur complement fieldsplit, determine how the Schur complement will be
2289: preconditioned. See PCFieldSplitSetSchurPre() for details.
2291: Logically Collective on PC
2293: Input Parameters:
2294: . pc - the preconditioner context
2296: Output Parameters:
2297: + ptype - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11, PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_PRE_USER
2298: - userpre - matrix to use for preconditioning (with PC_FIELDSPLIT_PRE_USER), or NULL
2300: Level: intermediate
2302: .seealso: PCFieldSplitSetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC
2304: @*/
2305: PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
2306: {
2311: PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre));
2312: return(0);
2313: }
2315: /*@
2316: PCFieldSplitSchurGetS - extract the MatSchurComplement object used by this PC in case it needs to be configured separately
2318: Not collective
2320: Input Parameter:
2321: . pc - the preconditioner context
2323: Output Parameter:
2324: . S - the Schur complement matrix
2326: Notes:
2327: This matrix should not be destroyed using MatDestroy(); rather, use PCFieldSplitSchurRestoreS().
2329: Level: advanced
2331: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurRestoreS()
2333: @*/
2334: PetscErrorCode PCFieldSplitSchurGetS(PC pc,Mat *S)
2335: {
2337: const char* t;
2338: PetscBool isfs;
2339: PC_FieldSplit *jac;
2343: PetscObjectGetType((PetscObject)pc,&t);
2344: PetscStrcmp(t,PCFIELDSPLIT,&isfs);
2345: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
2346: jac = (PC_FieldSplit*)pc->data;
2347: if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
2348: if (S) *S = jac->schur;
2349: return(0);
2350: }
2352: /*@
2353: PCFieldSplitSchurRestoreS - restores the MatSchurComplement object used by this PC
2355: Not collective
2357: Input Parameters:
2358: + pc - the preconditioner context
2359: . S - the Schur complement matrix
2361: Level: advanced
2363: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurGetS()
2365: @*/
2366: PetscErrorCode PCFieldSplitSchurRestoreS(PC pc,Mat *S)
2367: {
2369: const char* t;
2370: PetscBool isfs;
2371: PC_FieldSplit *jac;
2375: PetscObjectGetType((PetscObject)pc,&t);
2376: PetscStrcmp(t,PCFIELDSPLIT,&isfs);
2377: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
2378: jac = (PC_FieldSplit*)pc->data;
2379: if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
2380: if (!S || *S != jac->schur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten");
2381: return(0);
2382: }
2385: static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
2386: {
2387: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2391: jac->schurpre = ptype;
2392: if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) {
2393: MatDestroy(&jac->schur_user);
2394: jac->schur_user = pre;
2395: PetscObjectReference((PetscObject)jac->schur_user);
2396: }
2397: return(0);
2398: }
2400: static PetscErrorCode PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
2401: {
2402: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2405: *ptype = jac->schurpre;
2406: *pre = jac->schur_user;
2407: return(0);
2408: }
2410: /*@
2411: PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain in the preconditioner
2413: Collective on PC
2415: Input Parameters:
2416: + pc - the preconditioner context
2417: - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default
2419: Options Database:
2420: . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full
2423: Level: intermediate
2425: Notes:
2426: The FULL factorization is
2428: $ (A B) = (1 0) (A 0) (1 Ainv*B) = L D U
2429: $ (C E) (C*Ainv 1) (0 S) (0 1 )
2431: where S = E - C*Ainv*B. In practice, the full factorization is applied via block triangular solves with the grouping L*(D*U). UPPER uses D*U, LOWER uses L*D,
2432: and DIAG is the diagonal part with the sign of S flipped (because this makes the preconditioner positive definite for many formulations, thus allowing the use of KSPMINRES). Sign flipping of S can be turned off with PCFieldSplitSetSchurScale().
2434: $ If A and S are solved exactly
2435: $ *) FULL factorization is a direct solver.
2436: $ *) The preconditioned operator with LOWER or UPPER has all eigenvalues equal to 1 and minimal polynomial of degree 2, so KSPGMRES converges in 2 iterations.
2437: $ *) With DIAG, the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so KSPGMRES converges in at most 4 iterations.
2439: If the iteration count is very low, consider using KSPFGMRES or KSPGCR which can use one less preconditioner
2440: application in this case. Note that the preconditioned operator may be highly non-normal, so such fast convergence may not be observed in practice.
2442: For symmetric problems in which A is positive definite and S is negative definite, DIAG can be used with KSPMINRES.
2444: Note that a flexible method like KSPFGMRES or KSPGCR must be used if the fieldsplit preconditioner is nonlinear (e.g. a few iterations of a Krylov method is used to solve with A or S).
2446: References:
2447: + 1. - Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000).
2448: - 2. - Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001).
2450: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCFieldSplitSetSchurScale()
2451: @*/
2452: PetscErrorCode PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype)
2453: {
2458: PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));
2459: return(0);
2460: }
2462: static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype)
2463: {
2464: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2467: jac->schurfactorization = ftype;
2468: return(0);
2469: }
2471: /*@
2472: PCFieldSplitSetSchurScale - Controls the sign flip of S for PC_FIELDSPLIT_SCHUR_FACT_DIAG.
2474: Collective on PC
2476: Input Parameters:
2477: + pc - the preconditioner context
2478: - scale - scaling factor for the Schur complement
2480: Options Database:
2481: . -pc_fieldsplit_schur_scale - default is -1.0
2483: Level: intermediate
2485: .seealso: PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurFactType, PCFieldSplitSetSchurScale()
2486: @*/
2487: PetscErrorCode PCFieldSplitSetSchurScale(PC pc,PetscScalar scale)
2488: {
2494: PetscTryMethod(pc,"PCFieldSplitSetSchurScale_C",(PC,PetscScalar),(pc,scale));
2495: return(0);
2496: }
2498: static PetscErrorCode PCFieldSplitSetSchurScale_FieldSplit(PC pc,PetscScalar scale)
2499: {
2500: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2503: jac->schurscale = scale;
2504: return(0);
2505: }
2507: /*@C
2508: PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement
2510: Collective on KSP
2512: Input Parameter:
2513: . pc - the preconditioner context
2515: Output Parameters:
2516: + A00 - the (0,0) block
2517: . A01 - the (0,1) block
2518: . A10 - the (1,0) block
2519: - A11 - the (1,1) block
2521: Level: advanced
2523: .seealso: PCFIELDSPLIT
2524: @*/
2525: PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11)
2526: {
2527: PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;
2531: if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach.");
2532: if (A00) *A00 = jac->pmat[0];
2533: if (A01) *A01 = jac->B;
2534: if (A10) *A10 = jac->C;
2535: if (A11) *A11 = jac->pmat[1];
2536: return(0);
2537: }
2539: /*@
2540: PCFieldSplitSetGKBTol - Sets the solver tolerance for the generalized Golub-Kahan bidiagonalization preconditioner.
2542: Collective on PC
2544: Notes:
2545: The generalized GKB algorithm uses a lower bound estimate of the error in energy norm as stopping criterion.
2546: It stops once the lower bound estimate undershoots the required solver tolerance. Although the actual error might be bigger than
2547: this estimate, the stopping criterion is satisfactory in practical cases [A13].
2549: [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013.
2551: Input Parameters:
2552: + pc - the preconditioner context
2553: - tolerance - the solver tolerance
2555: Options Database:
2556: . -pc_fieldsplit_gkb_tol - default is 1e-5
2558: Level: intermediate
2560: .seealso: PCFIELDSPLIT, PCFieldSplitSetGKBDelay(), PCFieldSplitSetGKBNu(), PCFieldSplitSetGKBMaxit()
2561: @*/
2562: PetscErrorCode PCFieldSplitSetGKBTol(PC pc,PetscReal tolerance)
2563: {
2569: PetscTryMethod(pc,"PCFieldSplitSetGKBTol_C",(PC,PetscReal),(pc,tolerance));
2570: return(0);
2571: }
2573: static PetscErrorCode PCFieldSplitSetGKBTol_FieldSplit(PC pc,PetscReal tolerance)
2574: {
2575: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2578: jac->gkbtol = tolerance;
2579: return(0);
2580: }
2583: /*@
2584: PCFieldSplitSetGKBMaxit - Sets the maximum number of iterations for the generalized Golub-Kahan bidiagonalization
2585: preconditioner.
2587: Collective on PC
2589: Input Parameters:
2590: + pc - the preconditioner context
2591: - maxit - the maximum number of iterations
2593: Options Database:
2594: . -pc_fieldsplit_gkb_maxit - default is 100
2596: Level: intermediate
2598: .seealso: PCFIELDSPLIT, PCFieldSplitSetGKBDelay(), PCFieldSplitSetGKBTol(), PCFieldSplitSetGKBNu()
2599: @*/
2600: PetscErrorCode PCFieldSplitSetGKBMaxit(PC pc,PetscInt maxit)
2601: {
2607: PetscTryMethod(pc,"PCFieldSplitSetGKBMaxit_C",(PC,PetscInt),(pc,maxit));
2608: return(0);
2609: }
2611: static PetscErrorCode PCFieldSplitSetGKBMaxit_FieldSplit(PC pc,PetscInt maxit)
2612: {
2613: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2616: jac->gkbmaxit = maxit;
2617: return(0);
2618: }
2620: /*@
2621: PCFieldSplitSetGKBDelay - Sets the delay in the lower bound error estimate in the generalized Golub-Kahan bidiagonalization
2622: preconditioner.
2624: Collective on PC
2626: Notes:
2627: The algorithm uses a lower bound estimate of the error in energy norm as stopping criterion. The lower bound of the error ||u-u^k||_H
2628: is expressed as a truncated sum. The error at iteration k can only be measured at iteration (k + delay), and thus the algorithm needs
2629: at least (delay + 1) iterations to stop. For more details on the generalized Golub-Kahan bidiagonalization method and its lower bound stopping criterion, please refer to
2631: [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013.
2633: Input Parameters:
2634: + pc - the preconditioner context
2635: - delay - the delay window in the lower bound estimate
2637: Options Database:
2638: . -pc_fieldsplit_gkb_delay - default is 5
2640: Level: intermediate
2642: .seealso: PCFIELDSPLIT, PCFieldSplitSetGKBNu(), PCFieldSplitSetGKBTol(), PCFieldSplitSetGKBMaxit()
2643: @*/
2644: PetscErrorCode PCFieldSplitSetGKBDelay(PC pc,PetscInt delay)
2645: {
2651: PetscTryMethod(pc,"PCFieldSplitSetGKBDelay_C",(PC,PetscInt),(pc,delay));
2652: return(0);
2653: }
2655: static PetscErrorCode PCFieldSplitSetGKBDelay_FieldSplit(PC pc,PetscInt delay)
2656: {
2657: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2660: jac->gkbdelay = delay;
2661: return(0);
2662: }
2664: /*@
2665: PCFieldSplitSetGKBNu - Sets the scalar value nu >= 0 in the transformation H = A00 + nu*A01*A01' of the (1,1) block in the Golub-Kahan bidiagonalization preconditioner.
2667: Collective on PC
2669: Notes:
2670: This shift is in general done to obtain better convergence properties for the outer loop of the algorithm. This is often achieved by chosing nu sufficiently big. However,
2671: if nu is chosen too big, the matrix H might be badly conditioned and the solution of the linear system Hx = b in the inner loop gets difficult. It is therefore
2672: necessary to find a good balance in between the convergence of the inner and outer loop.
2674: For nu = 0, no shift is done. In this case A00 has to be positive definite. The matrix N in [Ar13] is then chosen as identity.
2676: [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013.
2678: Input Parameters:
2679: + pc - the preconditioner context
2680: - nu - the shift parameter
2682: Options Database:
2683: . -pc_fieldsplit_gkb_nu - default is 1
2685: Level: intermediate
2687: .seealso: PCFIELDSPLIT, PCFieldSplitSetGKBDelay(), PCFieldSplitSetGKBTol(), PCFieldSplitSetGKBMaxit()
2688: @*/
2689: PetscErrorCode PCFieldSplitSetGKBNu(PC pc,PetscReal nu)
2690: {
2696: PetscTryMethod(pc,"PCFieldSplitSetGKBNu_C",(PC,PetscReal),(pc,nu));
2697: return(0);
2698: }
2700: static PetscErrorCode PCFieldSplitSetGKBNu_FieldSplit(PC pc,PetscReal nu)
2701: {
2702: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2705: jac->gkbnu = nu;
2706: return(0);
2707: }
2710: static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type)
2711: {
2712: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2716: jac->type = type;
2718: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",0);
2719: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",0);
2720: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",0);
2721: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);
2722: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",0);
2723: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBTol_C",0);
2724: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBMaxit_C",0);
2725: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBNu_C",0);
2726: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBDelay_C",0);
2728: if (type == PC_COMPOSITE_SCHUR) {
2729: pc->ops->apply = PCApply_FieldSplit_Schur;
2730: pc->ops->view = PCView_FieldSplit_Schur;
2732: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);
2733: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit);
2734: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit);
2735: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);
2736: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",PCFieldSplitSetSchurScale_FieldSplit);
2737: } else if (type == PC_COMPOSITE_GKB){
2738: pc->ops->apply = PCApply_FieldSplit_GKB;
2739: pc->ops->view = PCView_FieldSplit_GKB;
2741: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2742: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBTol_C",PCFieldSplitSetGKBTol_FieldSplit);
2743: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBMaxit_C",PCFieldSplitSetGKBMaxit_FieldSplit);
2744: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBNu_C",PCFieldSplitSetGKBNu_FieldSplit);
2745: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBDelay_C",PCFieldSplitSetGKBDelay_FieldSplit);
2746: } else {
2747: pc->ops->apply = PCApply_FieldSplit;
2748: pc->ops->view = PCView_FieldSplit;
2750: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2751: }
2752: return(0);
2753: }
2755: static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs)
2756: {
2757: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2760: if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs);
2761: 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);
2762: jac->bs = bs;
2763: return(0);
2764: }
2766: /*@
2767: PCFieldSplitSetType - Sets the type of fieldsplit preconditioner.
2769: Collective on PC
2771: Input Parameter:
2772: . pc - the preconditioner context
2773: . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR
2775: Options Database Key:
2776: . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type
2778: Level: Intermediate
2780: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
2782: .seealso: PCCompositeSetType()
2784: @*/
2785: PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type)
2786: {
2791: PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));
2792: return(0);
2793: }
2795: /*@
2796: PCFieldSplitGetType - Gets the type of fieldsplit preconditioner.
2798: Not collective
2800: Input Parameter:
2801: . pc - the preconditioner context
2803: Output Parameter:
2804: . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR
2806: Level: Intermediate
2808: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
2809: .seealso: PCCompositeSetType()
2810: @*/
2811: PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type)
2812: {
2813: PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;
2818: *type = jac->type;
2819: return(0);
2820: }
2822: /*@
2823: PCFieldSplitSetDMSplits - Flags whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible.
2825: Logically Collective
2827: Input Parameters:
2828: + pc - the preconditioner context
2829: - flg - boolean indicating whether to use field splits defined by the DM
2831: Options Database Key:
2832: . -pc_fieldsplit_dm_splits
2834: Level: Intermediate
2836: .keywords: PC, DM, composite preconditioner, additive, multiplicative
2838: .seealso: PCFieldSplitGetDMSplits()
2840: @*/
2841: PetscErrorCode PCFieldSplitSetDMSplits(PC pc,PetscBool flg)
2842: {
2843: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2844: PetscBool isfs;
2850: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2851: if (isfs) {
2852: jac->dm_splits = flg;
2853: }
2854: return(0);
2855: }
2858: /*@
2859: PCFieldSplitGetDMSplits - Returns flag indicating whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible.
2861: Logically Collective
2863: Input Parameter:
2864: . pc - the preconditioner context
2866: Output Parameter:
2867: . flg - boolean indicating whether to use field splits defined by the DM
2869: Level: Intermediate
2871: .keywords: PC, DM, composite preconditioner, additive, multiplicative
2873: .seealso: PCFieldSplitSetDMSplits()
2875: @*/
2876: PetscErrorCode PCFieldSplitGetDMSplits(PC pc,PetscBool* flg)
2877: {
2878: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2879: PetscBool isfs;
2885: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2886: if (isfs) {
2887: if(flg) *flg = jac->dm_splits;
2888: }
2889: return(0);
2890: }
2892: /*@
2893: PCFieldSplitGetDetectSaddlePoint - Returns flag indicating whether PCFieldSplit will attempt to automatically determine fields based on zero diagonal entries.
2895: Logically Collective
2897: Input Parameter:
2898: . pc - the preconditioner context
2900: Output Parameter:
2901: . flg - boolean indicating whether to detect fields or not
2903: Level: Intermediate
2905: .seealso: PCFIELDSPLIT, PCFieldSplitSetDetectSaddlePoint()
2907: @*/
2908: PetscErrorCode PCFieldSplitGetDetectSaddlePoint(PC pc,PetscBool *flg)
2909: {
2910: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2913: *flg = jac->detect;
2914: return(0);
2915: }
2917: /*@
2918: PCFieldSplitSetDetectSaddlePoint - Sets flag indicating whether PCFieldSplit will attempt to automatically determine fields based on zero diagonal entries.
2920: Logically Collective
2922: Notes:
2923: Also sets the split type to PC_COMPOSITE_SCHUR (see PCFieldSplitSetType()) and the Schur preconditioner type to PC_FIELDSPLIT_SCHUR_PRE_SELF (see PCFieldSplitSetSchurPre()).
2925: Input Parameter:
2926: . pc - the preconditioner context
2928: Output Parameter:
2929: . flg - boolean indicating whether to detect fields or not
2931: Options Database Key:
2932: . -pc_fieldsplit_detect_saddle_point
2934: Level: Intermediate
2936: .seealso: PCFIELDSPLIT, PCFieldSplitSetDetectSaddlePoint(), PCFieldSplitSetType(), PCFieldSplitSetSchurPre()
2938: @*/
2939: PetscErrorCode PCFieldSplitSetDetectSaddlePoint(PC pc,PetscBool flg)
2940: {
2941: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2945: jac->detect = flg;
2946: if (jac->detect) {
2947: PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);
2948: PCFieldSplitSetSchurPre(pc,PC_FIELDSPLIT_SCHUR_PRE_SELF,NULL);
2949: }
2950: return(0);
2951: }
2953: /* -------------------------------------------------------------------------------------*/
2954: /*MC
2955: PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual
2956: fields or groups of fields. See the users manual section "Solving Block Matrices" for more details.
2958: To set options on the solvers for each block append -fieldsplit_ to all the PC
2959: options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1
2961: To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP()
2962: and set the options directly on the resulting KSP object
2964: Level: intermediate
2966: Options Database Keys:
2967: + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split
2968: . -pc_fieldsplit_default - automatically add any fields to additional splits that have not
2969: been supplied explicitly by -pc_fieldsplit_%d_fields
2970: . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields)
2971: . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur,gkb> - type of relaxation or factorization splitting
2972: . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11; see PCFieldSplitSetSchurPre()
2973: . -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero diagonal and uses Schur complement with no preconditioner as the solver
2975: . Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_
2976: for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields
2977: - Options prefix for inner solver when using Golub Kahan biadiagonalization preconditioner is -fieldsplit_0_
2979: Notes:
2980: Use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS()
2981: to define a field by an arbitrary collection of entries.
2983: If no fields are set the default is used. The fields are defined by entries strided by bs,
2984: beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(),
2985: if this is not called the block size defaults to the blocksize of the second matrix passed
2986: to KSPSetOperators()/PCSetOperators().
2988: $ For the Schur complement preconditioner if J = ( A00 A01 )
2989: $ ( A10 A11 )
2990: $ the preconditioner using full factorization is
2991: $ ( I -ksp(A00) A01 ) ( inv(A00) 0 ) ( I 0 )
2992: $ ( 0 I ) ( 0 ksp(S) ) ( -A10 ksp(A00) I )
2993: where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. S is the Schur complement
2994: $ S = A11 - A10 ksp(A00) A01
2995: 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
2996: in providing the SECOND split or 1 if not give). For PCFieldSplitGetSubKSP() when field number is 0,
2997: it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default
2998: A11 is used to construct a preconditioner for S, use PCFieldSplitSetSchurPre() for all the possible ways to construct the preconditioner for S.
3000: The factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above,
3001: diag gives
3002: $ ( inv(A00) 0 )
3003: $ ( 0 -ksp(S) )
3004: note that slightly counter intuitively there is a negative in front of the ksp(S) so that the preconditioner is positive definite. For SPD matrices J, the sign flip
3005: can be turned off with PCFieldSplitSetSchurScale() or by command line -pc_fieldsplit_schur_scale 1.0. The lower factorization is the inverse of
3006: $ ( A00 0 )
3007: $ ( A10 S )
3008: where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of
3009: $ ( A00 A01 )
3010: $ ( 0 S )
3011: where again the inverses of A00 and S are applied using KSPs.
3013: If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS
3014: is used automatically for a second block.
3016: The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1.
3017: Generally it should be used with the AIJ format.
3019: The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see,
3020: for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT
3021: inside a smoother resulting in "Distributive Smoothers".
3023: Concepts: physics based preconditioners, block preconditioners
3025: There is a nice discussion of block preconditioners in
3027: [El08] A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier-Stokes equations
3028: Howard Elman, V.E. Howle, John Shadid, Robert Shuttleworth, Ray Tuminaro, Journal of Computational Physics 227 (2008) 1790--1808
3029: http://chess.cs.umd.edu/~elman/papers/tax.pdf
3031: The Constrained Pressure Preconditioner (CPR) can be implemented using PCCOMPOSITE with PCGALERKIN. CPR first solves an R A P subsystem, updates the
3032: residual on all variables (PCCompositeSetType(pc,PC_COMPOSITE_MULTIPLICATIVE)), and then applies a simple ILU like preconditioner on all the variables.
3034: The generalized Golub-Kahan bidiagonalization preconditioner (gkb) can be applied to symmetric 2x2 block matrices of the shape
3035: $ ( A00 A01 )
3036: $ ( A01' 0 )
3037: with A00 positive semi-definite. The implementation follows [Ar13]. Therein, we choose N := 1/nu * I and the (1,1)-block of the matrix is modified to H = A00 + nu*A01*A01'.
3038: A linear system Hx = b has to be solved in each iteration of the GKB algorithm. This solver is chosen with the option prefix -fieldsplit_0_.
3040: [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013.
3042: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC,
3043: PCFieldSplitGetSubKSP(), PCFieldSplitSchurGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre(),
3044: MatSchurComplementSetAinvType(), PCFieldSplitSetSchurScale(),
3045: PCFieldSplitSetDetectSaddlePoint()
3046: M*/
3048: PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc)
3049: {
3051: PC_FieldSplit *jac;
3054: PetscNewLog(pc,&jac);
3056: jac->bs = -1;
3057: jac->nsplits = 0;
3058: jac->type = PC_COMPOSITE_MULTIPLICATIVE;
3059: jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */
3060: jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL;
3061: jac->schurscale = -1.0;
3062: jac->dm_splits = PETSC_TRUE;
3063: jac->detect = PETSC_FALSE;
3064: jac->gkbtol = 1e-5;
3065: jac->gkbdelay = 5;
3066: jac->gkbnu = 1;
3067: jac->gkbmaxit = 100;
3068: jac->gkbmonitor = PETSC_FALSE;
3070: pc->data = (void*)jac;
3072: pc->ops->apply = PCApply_FieldSplit;
3073: pc->ops->applytranspose = PCApplyTranspose_FieldSplit;
3074: pc->ops->setup = PCSetUp_FieldSplit;
3075: pc->ops->reset = PCReset_FieldSplit;
3076: pc->ops->destroy = PCDestroy_FieldSplit;
3077: pc->ops->setfromoptions = PCSetFromOptions_FieldSplit;
3078: pc->ops->view = PCView_FieldSplit;
3079: pc->ops->applyrichardson = 0;
3081: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurGetSubKSP_C",PCFieldSplitSchurGetSubKSP_FieldSplit);
3082: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
3083: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);
3084: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);
3085: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);
3086: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);
3087: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",PCFieldSplitRestrictIS_FieldSplit);
3088: return(0);
3089: }