Actual source code: fieldsplit.c
1: #include <petsc/private/pcimpl.h>
2: #include <petsc/private/kspimpl.h>
3: #include <petscdm.h>
5: const char *const PCFieldSplitSchurPreTypes[] = {"SELF","SELFP","A11","USER","FULL","PCFieldSplitSchurPreType","PC_FIELDSPLIT_SCHUR_PRE_",NULL};
6: const char *const PCFieldSplitSchurFactTypes[] = {"DIAG","LOWER","UPPER","FULL","PCFieldSplitSchurFactType","PC_FIELDSPLIT_SCHUR_FACT_",NULL};
8: 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;
10: typedef struct _PC_FieldSplitLink *PC_FieldSplitLink;
11: struct _PC_FieldSplitLink {
12: KSP ksp;
13: Vec x,y,z;
14: char *splitname;
15: PetscInt nfields;
16: PetscInt *fields,*fields_col;
17: VecScatter sctx;
18: IS is,is_col;
19: PC_FieldSplitLink next,previous;
20: PetscLogEvent event;
22: /* Used only when setting coordinates with PCSetCoordinates */
23: PetscInt dim;
24: PetscInt ndofs;
25: PetscReal *coords;
26: };
28: typedef struct {
29: PCCompositeType type;
30: PetscBool defaultsplit; /* Flag for a system with a set of 'k' scalar fields with the same layout (and bs = k) */
31: PetscBool splitdefined; /* Flag is set after the splits have been defined, to prevent more splits from being added */
32: PetscInt bs; /* Block size for IS and Mat structures */
33: PetscInt nsplits; /* Number of field divisions defined */
34: Vec *x,*y,w1,w2;
35: Mat *mat; /* The diagonal block for each split */
36: Mat *pmat; /* The preconditioning diagonal block for each split */
37: Mat *Afield; /* The rows of the matrix associated with each split */
38: PetscBool issetup;
40: /* Only used when Schur complement preconditioning is used */
41: Mat B; /* The (0,1) block */
42: Mat C; /* The (1,0) block */
43: Mat schur; /* The Schur complement S = A11 - A10 A00^{-1} A01, the KSP here, kspinner, is H_1 in [El08] */
44: Mat schurp; /* Assembled approximation to S built by MatSchurComplement to be used as a preconditioning matrix when solving with S */
45: Mat schur_user; /* User-provided preconditioning matrix for the Schur complement */
46: PCFieldSplitSchurPreType schurpre; /* Determines which preconditioning matrix is used for the Schur complement */
47: PCFieldSplitSchurFactType schurfactorization;
48: KSP kspschur; /* The solver for S */
49: KSP kspupper; /* The solver for A in the upper diagonal part of the factorization (H_2 in [El08]) */
50: PetscScalar schurscale; /* Scaling factor for the Schur complement solution with DIAG factorization */
52: /* Only used when Golub-Kahan bidiagonalization preconditioning is used */
53: Mat H; /* The modified matrix H = A00 + nu*A01*A01' */
54: PetscReal gkbtol; /* Stopping tolerance for lower bound estimate */
55: PetscInt gkbdelay; /* The delay window for the stopping criterion */
56: PetscReal gkbnu; /* Parameter for augmented Lagrangian H = A + nu*A01*A01' */
57: PetscInt gkbmaxit; /* Maximum number of iterations for outer loop */
58: PetscBool gkbmonitor; /* Monitor for gkb iterations and the lower bound error */
59: PetscViewer gkbviewer; /* Viewer context for gkbmonitor */
60: Vec u,v,d,Hu; /* Work vectors for the GKB algorithm */
61: PetscScalar *vecz; /* Contains intermediate values, eg for lower bound */
63: PC_FieldSplitLink head;
64: PetscBool isrestrict; /* indicates PCFieldSplitRestrictIS() has been last called on this object, hack */
65: PetscBool suboptionsset; /* Indicates that the KSPSetFromOptions() has been called on the sub-KSPs */
66: PetscBool dm_splits; /* Whether to use DMCreateFieldDecomposition() whenever possible */
67: PetscBool diag_use_amat; /* Whether to extract diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */
68: PetscBool offdiag_use_amat; /* Whether to extract off-diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */
69: PetscBool detect; /* Whether to form 2-way split by finding zero diagonal entries */
70: PetscBool coordinates_set; /* Whether PCSetCoordinates has been called */
71: } PC_FieldSplit;
73: /*
74: Notes:
75: there is no particular reason that pmat, x, and y are stored as arrays in PC_FieldSplit instead of
76: inside PC_FieldSplitLink, just historical. If you want to be able to add new fields after already using the
77: PC you could change this.
78: */
80: /* This helper is so that setting a user-provided preconditioning matrix is orthogonal to choosing to use it. This way the
81: * application-provided FormJacobian can provide this matrix without interfering with the user's (command-line) choices. */
82: static Mat FieldSplitSchurPre(PC_FieldSplit *jac)
83: {
84: switch (jac->schurpre) {
85: case PC_FIELDSPLIT_SCHUR_PRE_SELF: return jac->schur;
86: case PC_FIELDSPLIT_SCHUR_PRE_SELFP: return jac->schurp;
87: case PC_FIELDSPLIT_SCHUR_PRE_A11: return jac->pmat[1];
88: case PC_FIELDSPLIT_SCHUR_PRE_FULL: /* We calculate this and store it in schur_user */
89: case PC_FIELDSPLIT_SCHUR_PRE_USER: /* Use a user-provided matrix if it is given, otherwise diagonal block */
90: default:
91: return jac->schur_user ? jac->schur_user : jac->pmat[1];
92: }
93: }
95: #include <petscdraw.h>
96: static PetscErrorCode PCView_FieldSplit(PC pc,PetscViewer viewer)
97: {
98: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
99: PetscBool iascii,isdraw;
100: PetscInt i,j;
101: PC_FieldSplitLink ilink = jac->head;
103: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
104: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
105: if (iascii) {
106: if (jac->bs > 0) {
107: PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs);
108: } else {
109: PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D\n",PCCompositeTypes[jac->type],jac->nsplits);
110: }
111: if (pc->useAmat) {
112: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");
113: }
114: if (jac->diag_use_amat) {
115: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for diagonal blocks\n");
116: }
117: if (jac->offdiag_use_amat) {
118: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for off-diagonal blocks\n");
119: }
120: PetscViewerASCIIPrintf(viewer," Solver info for each split is in the following KSP objects:\n");
121: for (i=0; i<jac->nsplits; i++) {
122: if (ilink->fields) {
123: PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);
124: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
125: for (j=0; j<ilink->nfields; j++) {
126: if (j > 0) {
127: PetscViewerASCIIPrintf(viewer,",");
128: }
129: PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);
130: }
131: PetscViewerASCIIPrintf(viewer,"\n");
132: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
133: } else {
134: PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);
135: }
136: KSPView(ilink->ksp,viewer);
137: ilink = ilink->next;
138: }
139: }
141: if (isdraw) {
142: PetscDraw draw;
143: PetscReal x,y,w,wd;
145: PetscViewerDrawGetDraw(viewer,0,&draw);
146: PetscDrawGetCurrentPoint(draw,&x,&y);
147: w = 2*PetscMin(1.0 - x,x);
148: wd = w/(jac->nsplits + 1);
149: x = x - wd*(jac->nsplits-1)/2.0;
150: for (i=0; i<jac->nsplits; i++) {
151: PetscDrawPushCurrentPoint(draw,x,y);
152: KSPView(ilink->ksp,viewer);
153: PetscDrawPopCurrentPoint(draw);
154: x += wd;
155: ilink = ilink->next;
156: }
157: }
158: return 0;
159: }
161: static PetscErrorCode PCView_FieldSplit_Schur(PC pc,PetscViewer viewer)
162: {
163: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
164: PetscBool iascii,isdraw;
165: PetscInt i,j;
166: PC_FieldSplitLink ilink = jac->head;
167: MatSchurComplementAinvType atype;
169: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
170: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
171: if (iascii) {
172: if (jac->bs > 0) {
173: PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, blocksize = %D, factorization %s\n",jac->bs,PCFieldSplitSchurFactTypes[jac->schurfactorization]);
174: } else {
175: PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, factorization %s\n",PCFieldSplitSchurFactTypes[jac->schurfactorization]);
176: }
177: if (pc->useAmat) {
178: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");
179: }
180: switch (jac->schurpre) {
181: case PC_FIELDSPLIT_SCHUR_PRE_SELF:
182: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from S itself\n");
183: break;
184: case PC_FIELDSPLIT_SCHUR_PRE_SELFP:
185: MatSchurComplementGetAinvType(jac->schur,&atype);
186: 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;
187: case PC_FIELDSPLIT_SCHUR_PRE_A11:
188: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n");
189: break;
190: case PC_FIELDSPLIT_SCHUR_PRE_FULL:
191: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from the exact Schur complement\n");
192: break;
193: case PC_FIELDSPLIT_SCHUR_PRE_USER:
194: if (jac->schur_user) {
195: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from user provided matrix\n");
196: } else {
197: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n");
198: }
199: break;
200: default:
201: SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Invalid Schur preconditioning type: %d", jac->schurpre);
202: }
203: PetscViewerASCIIPrintf(viewer," Split info:\n");
204: PetscViewerASCIIPushTab(viewer);
205: for (i=0; i<jac->nsplits; i++) {
206: if (ilink->fields) {
207: PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);
208: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
209: for (j=0; j<ilink->nfields; j++) {
210: if (j > 0) {
211: PetscViewerASCIIPrintf(viewer,",");
212: }
213: PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);
214: }
215: PetscViewerASCIIPrintf(viewer,"\n");
216: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
217: } else {
218: PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);
219: }
220: ilink = ilink->next;
221: }
222: PetscViewerASCIIPrintf(viewer,"KSP solver for A00 block\n");
223: PetscViewerASCIIPushTab(viewer);
224: if (jac->head) {
225: KSPView(jac->head->ksp,viewer);
226: } else PetscViewerASCIIPrintf(viewer," not yet available\n");
227: PetscViewerASCIIPopTab(viewer);
228: if (jac->head && jac->kspupper != jac->head->ksp) {
229: PetscViewerASCIIPrintf(viewer,"KSP solver for upper A00 in upper triangular factor \n");
230: PetscViewerASCIIPushTab(viewer);
231: if (jac->kspupper) KSPView(jac->kspupper,viewer);
232: else PetscViewerASCIIPrintf(viewer," not yet available\n");
233: PetscViewerASCIIPopTab(viewer);
234: }
235: PetscViewerASCIIPrintf(viewer,"KSP solver for S = A11 - A10 inv(A00) A01 \n");
236: PetscViewerASCIIPushTab(viewer);
237: if (jac->kspschur) {
238: KSPView(jac->kspschur,viewer);
239: } else {
240: PetscViewerASCIIPrintf(viewer," not yet available\n");
241: }
242: PetscViewerASCIIPopTab(viewer);
243: PetscViewerASCIIPopTab(viewer);
244: } else if (isdraw && jac->head) {
245: PetscDraw draw;
246: PetscReal x,y,w,wd,h;
247: PetscInt cnt = 2;
248: char str[32];
250: PetscViewerDrawGetDraw(viewer,0,&draw);
251: PetscDrawGetCurrentPoint(draw,&x,&y);
252: if (jac->kspupper != jac->head->ksp) cnt++;
253: w = 2*PetscMin(1.0 - x,x);
254: wd = w/(cnt + 1);
256: PetscSNPrintf(str,32,"Schur fact. %s",PCFieldSplitSchurFactTypes[jac->schurfactorization]);
257: PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
258: y -= h;
259: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_USER && !jac->schur_user) {
260: PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[PC_FIELDSPLIT_SCHUR_PRE_A11]);
261: } else {
262: PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[jac->schurpre]);
263: }
264: PetscDrawStringBoxed(draw,x+wd*(cnt-1)/2.0,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
265: y -= h;
266: x = x - wd*(cnt-1)/2.0;
268: PetscDrawPushCurrentPoint(draw,x,y);
269: KSPView(jac->head->ksp,viewer);
270: PetscDrawPopCurrentPoint(draw);
271: if (jac->kspupper != jac->head->ksp) {
272: x += wd;
273: PetscDrawPushCurrentPoint(draw,x,y);
274: KSPView(jac->kspupper,viewer);
275: PetscDrawPopCurrentPoint(draw);
276: }
277: x += wd;
278: PetscDrawPushCurrentPoint(draw,x,y);
279: KSPView(jac->kspschur,viewer);
280: PetscDrawPopCurrentPoint(draw);
281: }
282: return 0;
283: }
285: static PetscErrorCode PCView_FieldSplit_GKB(PC pc,PetscViewer viewer)
286: {
287: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
288: PetscBool iascii,isdraw;
289: PetscInt i,j;
290: PC_FieldSplitLink ilink = jac->head;
292: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
293: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
294: if (iascii) {
295: if (jac->bs > 0) {
296: PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs);
297: } else {
298: PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D\n",PCCompositeTypes[jac->type],jac->nsplits);
299: }
300: if (pc->useAmat) {
301: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");
302: }
303: if (jac->diag_use_amat) {
304: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for diagonal blocks\n");
305: }
306: if (jac->offdiag_use_amat) {
307: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for off-diagonal blocks\n");
308: }
310: PetscViewerASCIIPrintf(viewer," Stopping tolerance=%.1e, delay in error estimate=%D, maximum iterations=%D\n",jac->gkbtol,jac->gkbdelay,jac->gkbmaxit);
311: PetscViewerASCIIPrintf(viewer," Solver info for H = A00 + nu*A01*A01' matrix:\n");
312: PetscViewerASCIIPushTab(viewer);
314: if (ilink->fields) {
315: PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",0);
316: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
317: for (j=0; j<ilink->nfields; j++) {
318: if (j > 0) {
319: PetscViewerASCIIPrintf(viewer,",");
320: }
321: PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);
322: }
323: PetscViewerASCIIPrintf(viewer,"\n");
324: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
325: } else {
326: PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",0);
327: }
328: KSPView(ilink->ksp,viewer);
330: PetscViewerASCIIPopTab(viewer);
331: }
333: if (isdraw) {
334: PetscDraw draw;
335: PetscReal x,y,w,wd;
337: PetscViewerDrawGetDraw(viewer,0,&draw);
338: PetscDrawGetCurrentPoint(draw,&x,&y);
339: w = 2*PetscMin(1.0 - x,x);
340: wd = w/(jac->nsplits + 1);
341: x = x - wd*(jac->nsplits-1)/2.0;
342: for (i=0; i<jac->nsplits; i++) {
343: PetscDrawPushCurrentPoint(draw,x,y);
344: KSPView(ilink->ksp,viewer);
345: PetscDrawPopCurrentPoint(draw);
346: x += wd;
347: ilink = ilink->next;
348: }
349: }
350: return 0;
351: }
353: /* Precondition: jac->bs is set to a meaningful value */
354: static PetscErrorCode PCFieldSplitSetRuntimeSplits_Private(PC pc)
355: {
356: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
357: PetscInt i,nfields,*ifields,nfields_col,*ifields_col;
358: PetscBool flg,flg_col;
359: char optionname[128],splitname[8],optionname_col[128];
361: PetscMalloc1(jac->bs,&ifields);
362: PetscMalloc1(jac->bs,&ifields_col);
363: for (i=0,flg=PETSC_TRUE;; i++) {
364: PetscSNPrintf(splitname,sizeof(splitname),"%D",i);
365: PetscSNPrintf(optionname,sizeof(optionname),"-pc_fieldsplit_%D_fields",i);
366: PetscSNPrintf(optionname_col,sizeof(optionname_col),"-pc_fieldsplit_%D_fields_col",i);
367: nfields = jac->bs;
368: nfields_col = jac->bs;
369: PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,optionname,ifields,&nfields,&flg);
370: PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,optionname_col,ifields_col,&nfields_col,&flg_col);
371: if (!flg) break;
372: else if (flg && !flg_col) {
374: PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields);
375: } else {
378: PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields_col);
379: }
380: }
381: if (i > 0) {
382: /* Makes command-line setting of splits take precedence over setting them in code.
383: Otherwise subsequent calls to PCFieldSplitSetIS() or PCFieldSplitSetFields() would
384: create new splits, which would probably not be what the user wanted. */
385: jac->splitdefined = PETSC_TRUE;
386: }
387: PetscFree(ifields);
388: PetscFree(ifields_col);
389: return 0;
390: }
392: static PetscErrorCode PCFieldSplitSetDefaults(PC pc)
393: {
394: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
395: PC_FieldSplitLink ilink = jac->head;
396: PetscBool fieldsplit_default = PETSC_FALSE,coupling = PETSC_FALSE;
397: PetscInt i;
399: /*
400: Kinda messy, but at least this now uses DMCreateFieldDecomposition().
401: Should probably be rewritten.
402: */
403: if (!ilink) {
404: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_coupling",&coupling,NULL);
405: if (pc->dm && jac->dm_splits && !jac->detect && !coupling) {
406: PetscInt numFields, f, i, j;
407: char **fieldNames;
408: IS *fields;
409: DM *dms;
410: DM subdm[128];
411: PetscBool flg;
413: DMCreateFieldDecomposition(pc->dm, &numFields, &fieldNames, &fields, &dms);
414: /* Allow the user to prescribe the splits */
415: for (i = 0, flg = PETSC_TRUE;; i++) {
416: PetscInt ifields[128];
417: IS compField;
418: char optionname[128], splitname[8];
419: PetscInt nfields = numFields;
421: PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%D_fields", i);
422: PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix, optionname, ifields, &nfields, &flg);
423: if (!flg) break;
425: DMCreateSubDM(pc->dm, nfields, ifields, &compField, &subdm[i]);
426: if (nfields == 1) {
427: PCFieldSplitSetIS(pc, fieldNames[ifields[0]], compField);
428: } else {
429: PetscSNPrintf(splitname, sizeof(splitname), "%D", i);
430: PCFieldSplitSetIS(pc, splitname, compField);
431: }
432: ISDestroy(&compField);
433: for (j = 0; j < nfields; ++j) {
434: f = ifields[j];
435: PetscFree(fieldNames[f]);
436: ISDestroy(&fields[f]);
437: }
438: }
439: if (i == 0) {
440: for (f = 0; f < numFields; ++f) {
441: PCFieldSplitSetIS(pc, fieldNames[f], fields[f]);
442: PetscFree(fieldNames[f]);
443: ISDestroy(&fields[f]);
444: }
445: } else {
446: for (j=0; j<numFields; j++) {
447: DMDestroy(dms+j);
448: }
449: PetscFree(dms);
450: PetscMalloc1(i, &dms);
451: for (j = 0; j < i; ++j) dms[j] = subdm[j];
452: }
453: PetscFree(fieldNames);
454: PetscFree(fields);
455: if (dms) {
456: PetscInfo(pc, "Setting up physics based fieldsplit preconditioner using the embedded DM\n");
457: for (ilink = jac->head, i = 0; ilink; ilink = ilink->next, ++i) {
458: const char *prefix;
459: PetscObjectGetOptionsPrefix((PetscObject)(ilink->ksp),&prefix);
460: PetscObjectSetOptionsPrefix((PetscObject)(dms[i]), prefix);
461: KSPSetDM(ilink->ksp, dms[i]);
462: KSPSetDMActive(ilink->ksp, PETSC_FALSE);
463: {
464: PetscErrorCode (*func)(KSP,Mat,Mat,void*);
465: void *ctx;
467: DMKSPGetComputeOperators(pc->dm, &func, &ctx);
468: DMKSPSetComputeOperators(dms[i], func, ctx);
469: }
470: PetscObjectIncrementTabLevel((PetscObject)dms[i],(PetscObject)ilink->ksp,0);
471: DMDestroy(&dms[i]);
472: }
473: PetscFree(dms);
474: }
475: } else {
476: if (jac->bs <= 0) {
477: if (pc->pmat) {
478: MatGetBlockSize(pc->pmat,&jac->bs);
479: } else jac->bs = 1;
480: }
482: if (jac->detect) {
483: IS zerodiags,rest;
484: PetscInt nmin,nmax;
486: MatGetOwnershipRange(pc->mat,&nmin,&nmax);
487: if (jac->diag_use_amat) {
488: MatFindZeroDiagonals(pc->mat,&zerodiags);
489: } else {
490: MatFindZeroDiagonals(pc->pmat,&zerodiags);
491: }
492: ISComplement(zerodiags,nmin,nmax,&rest);
493: PCFieldSplitSetIS(pc,"0",rest);
494: PCFieldSplitSetIS(pc,"1",zerodiags);
495: ISDestroy(&zerodiags);
496: ISDestroy(&rest);
497: } else if (coupling) {
498: IS coupling,rest;
499: PetscInt nmin,nmax;
501: MatGetOwnershipRange(pc->mat,&nmin,&nmax);
502: if (jac->offdiag_use_amat) {
503: MatFindOffBlockDiagonalEntries(pc->mat,&coupling);
504: } else {
505: MatFindOffBlockDiagonalEntries(pc->pmat,&coupling);
506: }
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: Mat M = pc->pmat;
523: PetscBool isnest;
525: PetscInfo(pc,"Using default splitting of fields\n");
526: PetscObjectTypeCompare((PetscObject)pc->pmat,MATNEST,&isnest);
527: if (!isnest) {
528: M = pc->mat;
529: PetscObjectTypeCompare((PetscObject)pc->mat,MATNEST,&isnest);
530: }
531: if (isnest) {
532: IS *fields;
533: PetscInt nf;
535: MatNestGetSize(M,&nf,NULL);
536: PetscMalloc1(nf,&fields);
537: MatNestGetISs(M,fields,NULL);
538: for (i=0;i<nf;i++) {
539: PCFieldSplitSetIS(pc,NULL,fields[i]);
540: }
541: PetscFree(fields);
542: } else {
543: for (i=0; i<jac->bs; i++) {
544: char splitname[8];
545: PetscSNPrintf(splitname,sizeof(splitname),"%D",i);
546: PCFieldSplitSetFields(pc,splitname,1,&i,&i);
547: }
548: jac->defaultsplit = PETSC_TRUE;
549: }
550: }
551: }
552: }
553: } else if (jac->nsplits == 1) {
554: if (ilink->is) {
555: IS is2;
556: PetscInt nmin,nmax;
558: MatGetOwnershipRange(pc->mat,&nmin,&nmax);
559: ISComplement(ilink->is,nmin,nmax,&is2);
560: PCFieldSplitSetIS(pc,"1",is2);
561: ISDestroy(&is2);
562: } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Must provide at least two sets of fields to PCFieldSplit()");
563: }
566: return 0;
567: }
569: static PetscErrorCode MatGolubKahanComputeExplicitOperator(Mat A,Mat B,Mat C,Mat *H,PetscReal gkbnu)
570: {
571: Mat BT,T;
572: PetscReal nrmT,nrmB;
574: MatHermitianTranspose(C,MAT_INITIAL_MATRIX,&T); /* Test if augmented matrix is symmetric */
575: MatAXPY(T,-1.0,B,DIFFERENT_NONZERO_PATTERN);
576: MatNorm(T,NORM_1,&nrmT);
577: MatNorm(B,NORM_1,&nrmB);
578: if (nrmB > 0) {
579: if (nrmT/nrmB >= PETSC_SMALL) {
580: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix is not symmetric/hermitian, GKB is not applicable.");
581: }
582: }
583: /* Compute augmented Lagrangian matrix H = A00 + nu*A01*A01'. This corresponds to */
584: /* setting N := 1/nu*I in [Ar13]. */
585: MatHermitianTranspose(B,MAT_INITIAL_MATRIX,&BT);
586: MatMatMult(B,BT,MAT_INITIAL_MATRIX,PETSC_DEFAULT,H); /* H = A01*A01' */
587: MatAYPX(*H,gkbnu,A,DIFFERENT_NONZERO_PATTERN); /* H = A00 + nu*A01*A01' */
589: MatDestroy(&BT);
590: MatDestroy(&T);
591: return 0;
592: }
594: PETSC_EXTERN PetscErrorCode PetscOptionsFindPairPrefix_Private(PetscOptions,const char pre[], const char name[],const char *value[],PetscBool *flg);
596: static PetscErrorCode PCSetUp_FieldSplit(PC pc)
597: {
598: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
599: PC_FieldSplitLink ilink;
600: PetscInt i,nsplit;
601: PetscBool sorted, sorted_col;
603: pc->failedreason = PC_NOERROR;
604: PCFieldSplitSetDefaults(pc);
605: nsplit = jac->nsplits;
606: ilink = jac->head;
608: /* get the matrices for each split */
609: if (!jac->issetup) {
610: PetscInt rstart,rend,nslots,bs;
612: jac->issetup = PETSC_TRUE;
614: /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */
615: if (jac->defaultsplit || !ilink->is) {
616: if (jac->bs <= 0) jac->bs = nsplit;
617: }
618: bs = jac->bs;
619: MatGetOwnershipRange(pc->pmat,&rstart,&rend);
620: nslots = (rend - rstart)/bs;
621: for (i=0; i<nsplit; i++) {
622: if (jac->defaultsplit) {
623: ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+i,nsplit,&ilink->is);
624: ISDuplicate(ilink->is,&ilink->is_col);
625: } else if (!ilink->is) {
626: if (ilink->nfields > 1) {
627: PetscInt *ii,*jj,j,k,nfields = ilink->nfields,*fields = ilink->fields,*fields_col = ilink->fields_col;
628: PetscMalloc1(ilink->nfields*nslots,&ii);
629: PetscMalloc1(ilink->nfields*nslots,&jj);
630: for (j=0; j<nslots; j++) {
631: for (k=0; k<nfields; k++) {
632: ii[nfields*j + k] = rstart + bs*j + fields[k];
633: jj[nfields*j + k] = rstart + bs*j + fields_col[k];
634: }
635: }
636: ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,ii,PETSC_OWN_POINTER,&ilink->is);
637: ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,jj,PETSC_OWN_POINTER,&ilink->is_col);
638: ISSetBlockSize(ilink->is, nfields);
639: ISSetBlockSize(ilink->is_col, nfields);
640: } else {
641: ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields[0],bs,&ilink->is);
642: ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields_col[0],bs,&ilink->is_col);
643: }
644: }
645: ISSorted(ilink->is,&sorted);
646: if (ilink->is_col) ISSorted(ilink->is_col,&sorted_col);
648: ilink = ilink->next;
649: }
650: }
652: ilink = jac->head;
653: if (!jac->pmat) {
654: Vec xtmp;
656: MatCreateVecs(pc->pmat,&xtmp,NULL);
657: PetscMalloc1(nsplit,&jac->pmat);
658: PetscMalloc2(nsplit,&jac->x,nsplit,&jac->y);
659: for (i=0; i<nsplit; i++) {
660: MatNullSpace sp;
662: /* Check for preconditioning matrix attached to IS */
663: PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &jac->pmat[i]);
664: if (jac->pmat[i]) {
665: PetscObjectReference((PetscObject) jac->pmat[i]);
666: if (jac->type == PC_COMPOSITE_SCHUR) {
667: jac->schur_user = jac->pmat[i];
669: PetscObjectReference((PetscObject) jac->schur_user);
670: }
671: } else {
672: const char *prefix;
673: MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->pmat[i]);
674: KSPGetOptionsPrefix(ilink->ksp,&prefix);
675: MatSetOptionsPrefix(jac->pmat[i],prefix);
676: MatViewFromOptions(jac->pmat[i],NULL,"-mat_view");
677: }
678: /* create work vectors for each split */
679: MatCreateVecs(jac->pmat[i],&jac->x[i],&jac->y[i]);
680: ilink->x = jac->x[i]; ilink->y = jac->y[i]; ilink->z = NULL;
681: /* compute scatter contexts needed by multiplicative versions and non-default splits */
682: VecScatterCreate(xtmp,ilink->is,jac->x[i],NULL,&ilink->sctx);
683: PetscObjectQuery((PetscObject) ilink->is, "nearnullspace", (PetscObject*) &sp);
684: if (sp) {
685: MatSetNearNullSpace(jac->pmat[i], sp);
686: }
687: ilink = ilink->next;
688: }
689: VecDestroy(&xtmp);
690: } else {
691: MatReuse scall;
692: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
693: for (i=0; i<nsplit; i++) {
694: MatDestroy(&jac->pmat[i]);
695: }
696: scall = MAT_INITIAL_MATRIX;
697: } else scall = MAT_REUSE_MATRIX;
699: for (i=0; i<nsplit; i++) {
700: Mat pmat;
702: /* Check for preconditioning matrix attached to IS */
703: PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &pmat);
704: if (!pmat) {
705: MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is_col,scall,&jac->pmat[i]);
706: }
707: ilink = ilink->next;
708: }
709: }
710: if (jac->diag_use_amat) {
711: ilink = jac->head;
712: if (!jac->mat) {
713: PetscMalloc1(nsplit,&jac->mat);
714: for (i=0; i<nsplit; i++) {
715: MatCreateSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->mat[i]);
716: ilink = ilink->next;
717: }
718: } else {
719: MatReuse scall;
720: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
721: for (i=0; i<nsplit; i++) {
722: MatDestroy(&jac->mat[i]);
723: }
724: scall = MAT_INITIAL_MATRIX;
725: } else scall = MAT_REUSE_MATRIX;
727: for (i=0; i<nsplit; i++) {
728: MatCreateSubMatrix(pc->mat,ilink->is,ilink->is_col,scall,&jac->mat[i]);
729: ilink = ilink->next;
730: }
731: }
732: } else {
733: jac->mat = jac->pmat;
734: }
736: /* Check for null space attached to IS */
737: ilink = jac->head;
738: for (i=0; i<nsplit; i++) {
739: MatNullSpace sp;
741: PetscObjectQuery((PetscObject) ilink->is, "nullspace", (PetscObject*) &sp);
742: if (sp) {
743: MatSetNullSpace(jac->mat[i], sp);
744: }
745: ilink = ilink->next;
746: }
748: if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR && jac->type != PC_COMPOSITE_GKB) {
749: /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */
750: /* FIXME: Can/should we reuse jac->mat whenever (jac->diag_use_amat) is true? */
751: ilink = jac->head;
752: if (nsplit == 2 && jac->type == PC_COMPOSITE_MULTIPLICATIVE) {
753: /* 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 */
754: if (!jac->Afield) {
755: PetscCalloc1(nsplit,&jac->Afield);
756: if (jac->offdiag_use_amat) {
757: MatCreateSubMatrix(pc->mat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);
758: } else {
759: MatCreateSubMatrix(pc->pmat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);
760: }
761: } else {
762: MatReuse scall;
764: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
765: MatDestroy(&jac->Afield[1]);
766: scall = MAT_INITIAL_MATRIX;
767: } else scall = MAT_REUSE_MATRIX;
769: if (jac->offdiag_use_amat) {
770: MatCreateSubMatrix(pc->mat,ilink->next->is,ilink->is,scall,&jac->Afield[1]);
771: } else {
772: MatCreateSubMatrix(pc->pmat,ilink->next->is,ilink->is,scall,&jac->Afield[1]);
773: }
774: }
775: } else {
776: if (!jac->Afield) {
777: PetscMalloc1(nsplit,&jac->Afield);
778: for (i=0; i<nsplit; i++) {
779: if (jac->offdiag_use_amat) {
780: MatCreateSubMatrix(pc->mat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);
781: } else {
782: MatCreateSubMatrix(pc->pmat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);
783: }
784: ilink = ilink->next;
785: }
786: } else {
787: MatReuse scall;
788: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
789: for (i=0; i<nsplit; i++) {
790: MatDestroy(&jac->Afield[i]);
791: }
792: scall = MAT_INITIAL_MATRIX;
793: } else scall = MAT_REUSE_MATRIX;
795: for (i=0; i<nsplit; i++) {
796: if (jac->offdiag_use_amat) {
797: MatCreateSubMatrix(pc->mat,ilink->is,NULL,scall,&jac->Afield[i]);
798: } else {
799: MatCreateSubMatrix(pc->pmat,ilink->is,NULL,scall,&jac->Afield[i]);
800: }
801: ilink = ilink->next;
802: }
803: }
804: }
805: }
807: if (jac->type == PC_COMPOSITE_SCHUR) {
808: IS ccis;
809: PetscBool isspd;
810: PetscInt rstart,rend;
811: char lscname[256];
812: PetscObject LSC_L;
816: /* If pc->mat is SPD, don't scale by -1 the Schur complement */
817: if (jac->schurscale == (PetscScalar)-1.0) {
818: MatGetOption(pc->pmat,MAT_SPD,&isspd);
819: jac->schurscale = (isspd == PETSC_TRUE) ? 1.0 : -1.0;
820: }
822: /* When extracting off-diagonal submatrices, we take complements from this range */
823: MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);
825: if (jac->schur) {
826: KSP kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper;
827: MatReuse scall;
829: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
830: scall = MAT_INITIAL_MATRIX;
831: MatDestroy(&jac->B);
832: MatDestroy(&jac->C);
833: } else scall = MAT_REUSE_MATRIX;
835: MatSchurComplementGetKSP(jac->schur, &kspInner);
836: ilink = jac->head;
837: ISComplement(ilink->is_col,rstart,rend,&ccis);
838: if (jac->offdiag_use_amat) {
839: MatCreateSubMatrix(pc->mat,ilink->is,ccis,scall,&jac->B);
840: } else {
841: MatCreateSubMatrix(pc->pmat,ilink->is,ccis,scall,&jac->B);
842: }
843: ISDestroy(&ccis);
844: ilink = ilink->next;
845: ISComplement(ilink->is_col,rstart,rend,&ccis);
846: if (jac->offdiag_use_amat) {
847: MatCreateSubMatrix(pc->mat,ilink->is,ccis,scall,&jac->C);
848: } else {
849: MatCreateSubMatrix(pc->pmat,ilink->is,ccis,scall,&jac->C);
850: }
851: ISDestroy(&ccis);
852: MatSchurComplementUpdateSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);
853: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
854: MatDestroy(&jac->schurp);
855: MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);
856: }
857: if (kspA != kspInner) {
858: KSPSetOperators(kspA,jac->mat[0],jac->pmat[0]);
859: }
860: if (kspUpper != kspA) {
861: KSPSetOperators(kspUpper,jac->mat[0],jac->pmat[0]);
862: }
863: KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));
864: } else {
865: const char *Dprefix;
866: char schurprefix[256], schurmatprefix[256];
867: char schurtestoption[256];
868: MatNullSpace sp;
869: PetscBool flg;
870: KSP kspt;
872: /* extract the A01 and A10 matrices */
873: ilink = jac->head;
874: ISComplement(ilink->is_col,rstart,rend,&ccis);
875: if (jac->offdiag_use_amat) {
876: MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
877: } else {
878: MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
879: }
880: ISDestroy(&ccis);
881: ilink = ilink->next;
882: ISComplement(ilink->is_col,rstart,rend,&ccis);
883: if (jac->offdiag_use_amat) {
884: MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
885: } else {
886: MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
887: }
888: ISDestroy(&ccis);
890: /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */
891: MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur);
892: MatSetType(jac->schur,MATSCHURCOMPLEMENT);
893: MatSchurComplementSetSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);
894: PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
895: MatSetOptionsPrefix(jac->schur,schurmatprefix);
896: MatSchurComplementGetKSP(jac->schur,&kspt);
897: KSPSetOptionsPrefix(kspt,schurmatprefix);
899: /* Note: this is not true in general */
900: MatGetNullSpace(jac->mat[1], &sp);
901: if (sp) {
902: MatSetNullSpace(jac->schur, sp);
903: }
905: PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);
906: PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
907: if (flg) {
908: DM dmInner;
909: KSP kspInner;
910: PC pcInner;
912: MatSchurComplementGetKSP(jac->schur, &kspInner);
913: KSPReset(kspInner);
914: KSPSetOperators(kspInner,jac->mat[0],jac->pmat[0]);
915: PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
916: /* Indent this deeper to emphasize the "inner" nature of this solver. */
917: PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject) pc, 2);
918: PetscObjectIncrementTabLevel((PetscObject)kspInner->pc, (PetscObject) pc, 2);
919: KSPSetOptionsPrefix(kspInner, schurprefix);
921: /* Set DM for new solver */
922: KSPGetDM(jac->head->ksp, &dmInner);
923: KSPSetDM(kspInner, dmInner);
924: KSPSetDMActive(kspInner, PETSC_FALSE);
926: /* Defaults to PCKSP as preconditioner */
927: KSPGetPC(kspInner, &pcInner);
928: PCSetType(pcInner, PCKSP);
929: PCKSPSetKSP(pcInner, jac->head->ksp);
930: } else {
931: /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or
932: * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact,
933: * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for
934: * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make
935: * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used
936: * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */
937: KSPSetType(jac->head->ksp,KSPGMRES);
938: MatSchurComplementSetKSP(jac->schur,jac->head->ksp);
939: }
940: KSPSetOperators(jac->head->ksp,jac->mat[0],jac->pmat[0]);
941: KSPSetFromOptions(jac->head->ksp);
942: MatSetFromOptions(jac->schur);
944: PetscObjectTypeCompare((PetscObject)jac->schur, MATSCHURCOMPLEMENT, &flg);
945: if (flg) { /* Need to do this otherwise PCSetUp_KSP will overwrite the amat of jac->head->ksp */
946: KSP kspInner;
947: PC pcInner;
949: MatSchurComplementGetKSP(jac->schur, &kspInner);
950: KSPGetPC(kspInner, &pcInner);
951: PetscObjectTypeCompare((PetscObject)pcInner, PCKSP, &flg);
952: if (flg) {
953: KSP ksp;
955: PCKSPGetKSP(pcInner, &ksp);
956: if (ksp == jac->head->ksp) {
957: PCSetUseAmat(pcInner, PETSC_TRUE);
958: }
959: }
960: }
961: PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);
962: PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
963: if (flg) {
964: DM dmInner;
966: PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
967: KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper);
968: KSPSetErrorIfNotConverged(jac->kspupper,pc->erroriffailure);
969: KSPSetOptionsPrefix(jac->kspupper, schurprefix);
970: PetscObjectIncrementTabLevel((PetscObject)jac->kspupper, (PetscObject) pc, 1);
971: PetscObjectIncrementTabLevel((PetscObject)jac->kspupper->pc, (PetscObject) pc, 1);
972: KSPGetDM(jac->head->ksp, &dmInner);
973: KSPSetDM(jac->kspupper, dmInner);
974: KSPSetDMActive(jac->kspupper, PETSC_FALSE);
975: KSPSetFromOptions(jac->kspupper);
976: KSPSetOperators(jac->kspupper,jac->mat[0],jac->pmat[0]);
977: VecDuplicate(jac->head->x, &jac->head->z);
978: } else {
979: jac->kspupper = jac->head->ksp;
980: PetscObjectReference((PetscObject) jac->head->ksp);
981: }
983: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
984: MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);
985: }
986: KSPCreate(PetscObjectComm((PetscObject)pc),&jac->kspschur);
987: KSPSetErrorIfNotConverged(jac->kspschur,pc->erroriffailure);
988: PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);
989: PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);
990: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) {
991: PC pcschur;
992: KSPGetPC(jac->kspschur,&pcschur);
993: PCSetType(pcschur,PCNONE);
994: /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */
995: } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) {
996: MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user);
997: }
998: KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));
999: KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix);
1000: KSPSetOptionsPrefix(jac->kspschur, Dprefix);
1001: /* propagate DM */
1002: {
1003: DM sdm;
1004: KSPGetDM(jac->head->next->ksp, &sdm);
1005: if (sdm) {
1006: KSPSetDM(jac->kspschur, sdm);
1007: KSPSetDMActive(jac->kspschur, PETSC_FALSE);
1008: }
1009: }
1010: /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
1011: /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */
1012: KSPSetFromOptions(jac->kspschur);
1013: }
1014: MatAssemblyBegin(jac->schur,MAT_FINAL_ASSEMBLY);
1015: MatAssemblyEnd(jac->schur,MAT_FINAL_ASSEMBLY);
1017: /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */
1018: PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_L",ilink->splitname);
1019: PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);
1020: if (!LSC_L) PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);
1021: if (LSC_L) PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L);
1022: PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_Lp",ilink->splitname);
1023: PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);
1024: if (!LSC_L) PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);
1025: if (LSC_L) PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L);
1026: } else if (jac->type == PC_COMPOSITE_GKB) {
1027: IS ccis;
1028: PetscInt rstart,rend;
1032: ilink = jac->head;
1034: /* When extracting off-diagonal submatrices, we take complements from this range */
1035: MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);
1037: ISComplement(ilink->is_col,rstart,rend,&ccis);
1038: if (jac->offdiag_use_amat) {
1039: MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
1040: } else {
1041: MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
1042: }
1043: ISDestroy(&ccis);
1044: /* Create work vectors for GKB algorithm */
1045: VecDuplicate(ilink->x,&jac->u);
1046: VecDuplicate(ilink->x,&jac->Hu);
1047: VecDuplicate(ilink->x,&jac->w2);
1048: ilink = ilink->next;
1049: ISComplement(ilink->is_col,rstart,rend,&ccis);
1050: if (jac->offdiag_use_amat) {
1051: MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
1052: } else {
1053: MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
1054: }
1055: ISDestroy(&ccis);
1056: /* Create work vectors for GKB algorithm */
1057: VecDuplicate(ilink->x,&jac->v);
1058: VecDuplicate(ilink->x,&jac->d);
1059: VecDuplicate(ilink->x,&jac->w1);
1060: MatGolubKahanComputeExplicitOperator(jac->mat[0],jac->B,jac->C,&jac->H,jac->gkbnu);
1061: PetscCalloc1(jac->gkbdelay,&jac->vecz);
1063: ilink = jac->head;
1064: KSPSetOperators(ilink->ksp,jac->H,jac->H);
1065: if (!jac->suboptionsset) KSPSetFromOptions(ilink->ksp);
1066: /* Create gkb_monitor context */
1067: if (jac->gkbmonitor) {
1068: PetscInt tablevel;
1069: PetscViewerCreate(PETSC_COMM_WORLD,&jac->gkbviewer);
1070: PetscViewerSetType(jac->gkbviewer,PETSCVIEWERASCII);
1071: PetscObjectGetTabLevel((PetscObject)ilink->ksp,&tablevel);
1072: PetscViewerASCIISetTab(jac->gkbviewer,tablevel);
1073: PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)ilink->ksp,1);
1074: }
1075: } else {
1076: /* set up the individual splits' PCs */
1077: i = 0;
1078: ilink = jac->head;
1079: while (ilink) {
1080: KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i]);
1081: /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
1082: if (!jac->suboptionsset) KSPSetFromOptions(ilink->ksp);
1083: i++;
1084: ilink = ilink->next;
1085: }
1086: }
1088: /* Set coordinates to the sub PC objects whenever these are set */
1089: if (jac->coordinates_set) {
1090: PC pc_coords;
1091: if (jac->type == PC_COMPOSITE_SCHUR) {
1092: // Head is first block.
1093: KSPGetPC(jac->head->ksp, &pc_coords);
1094: PCSetCoordinates(pc_coords, jac->head->dim, jac->head->ndofs, jac->head->coords);
1095: // Second one is Schur block, but its KSP object is in kspschur.
1096: KSPGetPC(jac->kspschur, &pc_coords);
1097: PCSetCoordinates(pc_coords, jac->head->next->dim, jac->head->next->ndofs, jac->head->next->coords);
1098: } else if (jac->type == PC_COMPOSITE_GKB) {
1099: PetscInfo(pc, "Warning: Setting coordinates does nothing for the GKB Fieldpslit preconditioner");
1100: } else {
1101: ilink = jac->head;
1102: while (ilink) {
1103: KSPGetPC(ilink->ksp, &pc_coords);
1104: PCSetCoordinates(pc_coords, ilink->dim, ilink->ndofs, ilink->coords);
1105: ilink = ilink->next;
1106: }
1107: }
1108: }
1110: jac->suboptionsset = PETSC_TRUE;
1111: return 0;
1112: }
1114: #define FieldSplitSplitSolveAdd(ilink,xx,yy) \
1115: (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
1116: VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
1117: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\
1118: KSPSolve(ilink->ksp,ilink->x,ilink->y) || \
1119: KSPCheckSolve(ilink->ksp,pc,ilink->y) || \
1120: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\
1121: VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \
1122: VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE))
1124: static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y)
1125: {
1126: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1127: PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
1128: KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;
1130: switch (jac->schurfactorization) {
1131: case PC_FIELDSPLIT_SCHUR_FACT_DIAG:
1132: /* [A00 0; 0 -S], positive definite, suitable for MINRES */
1133: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1134: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
1135: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1136: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1137: KSPSolve(kspA,ilinkA->x,ilinkA->y);
1138: KSPCheckSolve(kspA,pc,ilinkA->y);
1139: PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1140: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1141: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
1142: PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1143: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
1144: KSPCheckSolve(jac->kspschur,pc,ilinkD->y);
1145: PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1146: VecScale(ilinkD->y,jac->schurscale);
1147: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1148: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1149: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1150: break;
1151: case PC_FIELDSPLIT_SCHUR_FACT_LOWER:
1152: /* [A00 0; A10 S], suitable for left preconditioning */
1153: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1154: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1155: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1156: KSPSolve(kspA,ilinkA->x,ilinkA->y);
1157: KSPCheckSolve(kspA,pc,ilinkA->y);
1158: PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1159: MatMult(jac->C,ilinkA->y,ilinkD->x);
1160: VecScale(ilinkD->x,-1.);
1161: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
1162: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1163: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
1164: PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1165: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
1166: KSPCheckSolve(jac->kspschur,pc,ilinkD->y);
1167: PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1168: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1169: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1170: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1171: break;
1172: case PC_FIELDSPLIT_SCHUR_FACT_UPPER:
1173: /* [A00 A01; 0 S], suitable for right preconditioning */
1174: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
1175: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
1176: PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1177: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
1178: KSPCheckSolve(jac->kspschur,pc,ilinkD->y);
1179: PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL)); PetscCall(MatMult(jac->B,ilinkD->y,ilinkA->x);
1180: VecScale(ilinkA->x,-1.);
1181: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);
1182: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1183: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);
1184: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1185: KSPSolve(kspA,ilinkA->x,ilinkA->y);
1186: KSPCheckSolve(kspA,pc,ilinkA->y);
1187: PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1188: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1189: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1190: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1191: break;
1192: case PC_FIELDSPLIT_SCHUR_FACT_FULL:
1193: /* [1 0; A10 A00^{-1} 1] [A00 0; 0 S] [1 A00^{-1}A01; 0 1] */
1194: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1195: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1196: PetscLogEventBegin(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);
1197: KSPSolve(kspLower,ilinkA->x,ilinkA->y);
1198: KSPCheckSolve(kspLower,pc,ilinkA->y);
1199: PetscLogEventEnd(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);
1200: MatMult(jac->C,ilinkA->y,ilinkD->x);
1201: VecScale(ilinkD->x,-1.0);
1202: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
1203: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
1205: PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1206: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
1207: KSPCheckSolve(jac->kspschur,pc,ilinkD->y);
1208: PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1209: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1211: if (kspUpper == kspA) {
1212: MatMult(jac->B,ilinkD->y,ilinkA->y);
1213: VecAXPY(ilinkA->x,-1.0,ilinkA->y);
1214: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1215: KSPSolve(kspA,ilinkA->x,ilinkA->y);
1216: KSPCheckSolve(kspA,pc,ilinkA->y);
1217: PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1218: } else {
1219: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1220: KSPSolve(kspA,ilinkA->x,ilinkA->y);
1221: KSPCheckSolve(kspA,pc,ilinkA->y);
1222: MatMult(jac->B,ilinkD->y,ilinkA->x);
1223: PetscLogEventBegin(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);
1224: KSPSolve(kspUpper,ilinkA->x,ilinkA->z);
1225: KSPCheckSolve(kspUpper,pc,ilinkA->z);
1226: PetscLogEventEnd(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);
1227: VecAXPY(ilinkA->y,-1.0,ilinkA->z);
1228: }
1229: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1230: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1231: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1232: }
1233: return 0;
1234: }
1236: static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y)
1237: {
1238: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1239: PC_FieldSplitLink ilink = jac->head;
1240: PetscInt cnt,bs;
1242: if (jac->type == PC_COMPOSITE_ADDITIVE) {
1243: if (jac->defaultsplit) {
1244: VecGetBlockSize(x,&bs);
1246: VecGetBlockSize(y,&bs);
1248: VecStrideGatherAll(x,jac->x,INSERT_VALUES);
1249: while (ilink) {
1250: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1251: KSPSolve(ilink->ksp,ilink->x,ilink->y);
1252: KSPCheckSolve(ilink->ksp,pc,ilink->y);
1253: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1254: ilink = ilink->next;
1255: }
1256: VecStrideScatterAll(jac->y,y,INSERT_VALUES);
1257: } else {
1258: VecSet(y,0.0);
1259: while (ilink) {
1260: FieldSplitSplitSolveAdd(ilink,x,y);
1261: ilink = ilink->next;
1262: }
1263: }
1264: } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) {
1265: VecSet(y,0.0);
1266: /* solve on first block for first block variables */
1267: VecScatterBegin(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);
1268: VecScatterEnd(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);
1269: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1270: KSPSolve(ilink->ksp,ilink->x,ilink->y);
1271: KSPCheckSolve(ilink->ksp,pc,ilink->y);
1272: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1273: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1274: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1276: /* compute the residual only onto second block variables using first block variables */
1277: MatMult(jac->Afield[1],ilink->y,ilink->next->x);
1278: ilink = ilink->next;
1279: VecScale(ilink->x,-1.0);
1280: VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1281: VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1283: /* solve on second block variables */
1284: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1285: KSPSolve(ilink->ksp,ilink->x,ilink->y);
1286: KSPCheckSolve(ilink->ksp,pc,ilink->y);
1287: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1288: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1289: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1290: } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1291: if (!jac->w1) {
1292: VecDuplicate(x,&jac->w1);
1293: VecDuplicate(x,&jac->w2);
1294: }
1295: VecSet(y,0.0);
1296: FieldSplitSplitSolveAdd(ilink,x,y);
1297: cnt = 1;
1298: while (ilink->next) {
1299: ilink = ilink->next;
1300: /* compute the residual only over the part of the vector needed */
1301: MatMult(jac->Afield[cnt++],y,ilink->x);
1302: VecScale(ilink->x,-1.0);
1303: VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1304: VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1305: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1306: KSPSolve(ilink->ksp,ilink->x,ilink->y);
1307: KSPCheckSolve(ilink->ksp,pc,ilink->y);
1308: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1309: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1310: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1311: }
1312: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1313: cnt -= 2;
1314: while (ilink->previous) {
1315: ilink = ilink->previous;
1316: /* compute the residual only over the part of the vector needed */
1317: MatMult(jac->Afield[cnt--],y,ilink->x);
1318: VecScale(ilink->x,-1.0);
1319: VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1320: VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1321: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1322: KSPSolve(ilink->ksp,ilink->x,ilink->y);
1323: KSPCheckSolve(ilink->ksp,pc,ilink->y);
1324: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1325: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1326: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1327: }
1328: }
1329: } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type);
1330: return 0;
1331: }
1333: static PetscErrorCode PCApply_FieldSplit_GKB(PC pc,Vec x,Vec y)
1334: {
1335: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1336: PC_FieldSplitLink ilinkA = jac->head,ilinkD = ilinkA->next;
1337: KSP ksp = ilinkA->ksp;
1338: Vec u,v,Hu,d,work1,work2;
1339: PetscScalar alpha,z,nrmz2,*vecz;
1340: PetscReal lowbnd,nu,beta;
1341: PetscInt j,iterGKB;
1343: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1344: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
1345: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1346: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
1348: u = jac->u;
1349: v = jac->v;
1350: Hu = jac->Hu;
1351: d = jac->d;
1352: work1 = jac->w1;
1353: work2 = jac->w2;
1354: vecz = jac->vecz;
1356: /* Change RHS to comply with matrix regularization H = A + nu*B*B' */
1357: /* Add q = q + nu*B*b */
1358: if (jac->gkbnu) {
1359: nu = jac->gkbnu;
1360: VecScale(ilinkD->x,jac->gkbnu);
1361: MatMultAdd(jac->B,ilinkD->x,ilinkA->x,ilinkA->x); /* q = q + nu*B*b */
1362: } else {
1363: /* Situation when no augmented Lagrangian is used. Then we set inner */
1364: /* matrix N = I in [Ar13], and thus nu = 1. */
1365: nu = 1;
1366: }
1368: /* Transform rhs from [q,tilde{b}] to [0,b] */
1369: PetscLogEventBegin(ilinkA->event,ksp,ilinkA->x,ilinkA->y,NULL);
1370: KSPSolve(ksp,ilinkA->x,ilinkA->y);
1371: KSPCheckSolve(ksp,pc,ilinkA->y);
1372: PetscLogEventEnd(ilinkA->event,ksp,ilinkA->x,ilinkA->y,NULL);
1373: MatMultHermitianTranspose(jac->B,ilinkA->y,work1);
1374: VecAXPBY(work1,1.0/nu,-1.0,ilinkD->x); /* c = b - B'*x */
1376: /* First step of algorithm */
1377: VecNorm(work1,NORM_2,&beta); /* beta = sqrt(nu*c'*c)*/
1378: KSPCheckDot(ksp,beta);
1379: beta = PetscSqrtReal(nu)*beta;
1380: VecAXPBY(v,nu/beta,0.0,work1); /* v = nu/beta *c */
1381: MatMult(jac->B,v,work2); /* u = H^{-1}*B*v */
1382: PetscLogEventBegin(ilinkA->event,ksp,work2,u,NULL);
1383: KSPSolve(ksp,work2,u);
1384: KSPCheckSolve(ksp,pc,u);
1385: PetscLogEventEnd(ilinkA->event,ksp,work2,u,NULL);
1386: MatMult(jac->H,u,Hu); /* alpha = u'*H*u */
1387: VecDot(Hu,u,&alpha);
1388: KSPCheckDot(ksp,alpha);
1390: alpha = PetscSqrtReal(PetscAbsScalar(alpha));
1391: VecScale(u,1.0/alpha);
1392: VecAXPBY(d,1.0/alpha,0.0,v); /* v = nu/beta *c */
1394: z = beta/alpha;
1395: vecz[1] = z;
1397: /* Computation of first iterate x(1) and p(1) */
1398: VecAXPY(ilinkA->y,z,u);
1399: VecCopy(d,ilinkD->y);
1400: VecScale(ilinkD->y,-z);
1402: iterGKB = 1; lowbnd = 2*jac->gkbtol;
1403: if (jac->gkbmonitor) {
1404: PetscViewerASCIIPrintf(jac->gkbviewer,"%3D GKB Lower bound estimate %14.12e\n",iterGKB,lowbnd);
1405: }
1407: while (iterGKB < jac->gkbmaxit && lowbnd > jac->gkbtol) {
1408: iterGKB += 1;
1409: MatMultHermitianTranspose(jac->B,u,work1); /* v <- nu*(B'*u-alpha/nu*v) */
1410: VecAXPBY(v,nu,-alpha,work1);
1411: VecNorm(v,NORM_2,&beta); /* beta = sqrt(nu)*v'*v */
1412: beta = beta/PetscSqrtReal(nu);
1413: VecScale(v,1.0/beta);
1414: MatMult(jac->B,v,work2); /* u <- H^{-1}*(B*v-beta*H*u) */
1415: MatMult(jac->H,u,Hu);
1416: VecAXPY(work2,-beta,Hu);
1417: PetscLogEventBegin(ilinkA->event,ksp,work2,u,NULL);
1418: KSPSolve(ksp,work2,u);
1419: KSPCheckSolve(ksp,pc,u);
1420: PetscLogEventEnd(ilinkA->event,ksp,work2,u,NULL);
1421: MatMult(jac->H,u,Hu); /* alpha = u'*H*u */
1422: VecDot(Hu,u,&alpha);
1423: KSPCheckDot(ksp,alpha);
1425: alpha = PetscSqrtReal(PetscAbsScalar(alpha));
1426: VecScale(u,1.0/alpha);
1428: z = -beta/alpha*z; /* z <- beta/alpha*z */
1429: vecz[0] = z;
1431: /* Computation of new iterate x(i+1) and p(i+1) */
1432: VecAXPBY(d,1.0/alpha,-beta/alpha,v); /* d = (v-beta*d)/alpha */
1433: VecAXPY(ilinkA->y,z,u); /* r = r + z*u */
1434: VecAXPY(ilinkD->y,-z,d); /* p = p - z*d */
1435: MatMult(jac->H,ilinkA->y,Hu); /* ||u||_H = u'*H*u */
1436: VecDot(Hu,ilinkA->y,&nrmz2);
1438: /* Compute Lower Bound estimate */
1439: if (iterGKB > jac->gkbdelay) {
1440: lowbnd = 0.0;
1441: for (j=0; j<jac->gkbdelay; j++) {
1442: lowbnd += PetscAbsScalar(vecz[j]*vecz[j]);
1443: }
1444: lowbnd = PetscSqrtReal(lowbnd/PetscAbsScalar(nrmz2));
1445: }
1447: for (j=0; j<jac->gkbdelay-1; j++) {
1448: vecz[jac->gkbdelay-j-1] = vecz[jac->gkbdelay-j-2];
1449: }
1450: if (jac->gkbmonitor) {
1451: PetscViewerASCIIPrintf(jac->gkbviewer,"%3D GKB Lower bound estimate %14.12e\n",iterGKB,lowbnd);
1452: }
1453: }
1455: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1456: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1457: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1458: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1460: return 0;
1461: }
1463: #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \
1464: (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
1465: VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
1466: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \
1467: KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \
1468: KSPCheckSolve(ilink->ksp,pc,ilink->x) || \
1469: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \
1470: VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \
1471: VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE))
1473: static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y)
1474: {
1475: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1476: PC_FieldSplitLink ilink = jac->head;
1477: PetscInt bs;
1479: if (jac->type == PC_COMPOSITE_ADDITIVE) {
1480: if (jac->defaultsplit) {
1481: VecGetBlockSize(x,&bs);
1483: VecGetBlockSize(y,&bs);
1485: VecStrideGatherAll(x,jac->x,INSERT_VALUES);
1486: while (ilink) {
1487: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1488: KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);
1489: KSPCheckSolve(ilink->ksp,pc,ilink->y);
1490: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1491: ilink = ilink->next;
1492: }
1493: VecStrideScatterAll(jac->y,y,INSERT_VALUES);
1494: } else {
1495: VecSet(y,0.0);
1496: while (ilink) {
1497: FieldSplitSplitSolveAddTranspose(ilink,x,y);
1498: ilink = ilink->next;
1499: }
1500: }
1501: } else {
1502: if (!jac->w1) {
1503: VecDuplicate(x,&jac->w1);
1504: VecDuplicate(x,&jac->w2);
1505: }
1506: VecSet(y,0.0);
1507: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1508: FieldSplitSplitSolveAddTranspose(ilink,x,y);
1509: while (ilink->next) {
1510: ilink = ilink->next;
1511: MatMultTranspose(pc->mat,y,jac->w1);
1512: VecWAXPY(jac->w2,-1.0,jac->w1,x);
1513: FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1514: }
1515: while (ilink->previous) {
1516: ilink = ilink->previous;
1517: MatMultTranspose(pc->mat,y,jac->w1);
1518: VecWAXPY(jac->w2,-1.0,jac->w1,x);
1519: FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1520: }
1521: } else {
1522: while (ilink->next) { /* get to last entry in linked list */
1523: ilink = ilink->next;
1524: }
1525: FieldSplitSplitSolveAddTranspose(ilink,x,y);
1526: while (ilink->previous) {
1527: ilink = ilink->previous;
1528: MatMultTranspose(pc->mat,y,jac->w1);
1529: VecWAXPY(jac->w2,-1.0,jac->w1,x);
1530: FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1531: }
1532: }
1533: }
1534: return 0;
1535: }
1537: static PetscErrorCode PCReset_FieldSplit(PC pc)
1538: {
1539: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1540: PC_FieldSplitLink ilink = jac->head,next;
1542: while (ilink) {
1543: KSPDestroy(&ilink->ksp);
1544: VecDestroy(&ilink->x);
1545: VecDestroy(&ilink->y);
1546: VecDestroy(&ilink->z);
1547: VecScatterDestroy(&ilink->sctx);
1548: ISDestroy(&ilink->is);
1549: ISDestroy(&ilink->is_col);
1550: PetscFree(ilink->splitname);
1551: PetscFree(ilink->fields);
1552: PetscFree(ilink->fields_col);
1553: next = ilink->next;
1554: PetscFree(ilink);
1555: ilink = next;
1556: }
1557: jac->head = NULL;
1558: PetscFree2(jac->x,jac->y);
1559: if (jac->mat && jac->mat != jac->pmat) {
1560: MatDestroyMatrices(jac->nsplits,&jac->mat);
1561: } else if (jac->mat) {
1562: jac->mat = NULL;
1563: }
1564: if (jac->pmat) MatDestroyMatrices(jac->nsplits,&jac->pmat);
1565: if (jac->Afield) MatDestroyMatrices(jac->nsplits,&jac->Afield);
1566: jac->nsplits = 0;
1567: VecDestroy(&jac->w1);
1568: VecDestroy(&jac->w2);
1569: MatDestroy(&jac->schur);
1570: MatDestroy(&jac->schurp);
1571: MatDestroy(&jac->schur_user);
1572: KSPDestroy(&jac->kspschur);
1573: KSPDestroy(&jac->kspupper);
1574: MatDestroy(&jac->B);
1575: MatDestroy(&jac->C);
1576: MatDestroy(&jac->H);
1577: VecDestroy(&jac->u);
1578: VecDestroy(&jac->v);
1579: VecDestroy(&jac->Hu);
1580: VecDestroy(&jac->d);
1581: PetscFree(jac->vecz);
1582: PetscViewerDestroy(&jac->gkbviewer);
1583: jac->isrestrict = PETSC_FALSE;
1584: return 0;
1585: }
1587: static PetscErrorCode PCDestroy_FieldSplit(PC pc)
1588: {
1589: PCReset_FieldSplit(pc);
1590: PetscFree(pc->data);
1591: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurGetSubKSP_C",NULL);
1592: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL);
1593: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL);
1594: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL);
1595: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL);
1596: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL);
1597: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL);
1598: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL);
1599: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL);
1600: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",NULL);
1601: return 0;
1602: }
1604: static PetscErrorCode PCSetFromOptions_FieldSplit(PetscOptionItems *PetscOptionsObject,PC pc)
1605: {
1606: PetscInt bs;
1607: PetscBool flg;
1608: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1609: PCCompositeType ctype;
1611: PetscOptionsHead(PetscOptionsObject,"FieldSplit options");
1612: PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL);
1613: PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);
1614: if (flg) {
1615: PCFieldSplitSetBlockSize(pc,bs);
1616: }
1617: jac->diag_use_amat = pc->useAmat;
1618: 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);
1619: jac->offdiag_use_amat = pc->useAmat;
1620: 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);
1621: PetscOptionsBool("-pc_fieldsplit_detect_saddle_point","Form 2-way split by detecting zero diagonal entries", "PCFieldSplitSetDetectSaddlePoint",jac->detect,&jac->detect,NULL);
1622: PCFieldSplitSetDetectSaddlePoint(pc,jac->detect); /* Sets split type and Schur PC type */
1623: PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);
1624: if (flg) {
1625: PCFieldSplitSetType(pc,ctype);
1626: }
1627: /* Only setup fields once */
1628: if ((jac->bs > 0) && (jac->nsplits == 0)) {
1629: /* only allow user to set fields from command line if bs is already known.
1630: otherwise user can set them in PCFieldSplitSetDefaults() */
1631: PCFieldSplitSetRuntimeSplits_Private(pc);
1632: if (jac->splitdefined) PetscInfo(pc,"Splits defined using the options database\n");
1633: }
1634: if (jac->type == PC_COMPOSITE_SCHUR) {
1635: PetscOptionsGetEnum(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);
1636: if (flg) PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");
1637: 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);
1638: PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSetSchurPre",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL);
1639: PetscOptionsScalar("-pc_fieldsplit_schur_scale","Scale Schur complement","PCFieldSplitSetSchurScale",jac->schurscale,&jac->schurscale,NULL);
1640: } else if (jac->type == PC_COMPOSITE_GKB) {
1641: PetscOptionsReal("-pc_fieldsplit_gkb_tol","The tolerance for the lower bound stopping criterion","PCFieldSplitGKBTol",jac->gkbtol,&jac->gkbtol,NULL);
1642: PetscOptionsInt("-pc_fieldsplit_gkb_delay","The delay value for lower bound criterion","PCFieldSplitGKBDelay",jac->gkbdelay,&jac->gkbdelay,NULL);
1643: PetscOptionsReal("-pc_fieldsplit_gkb_nu","Parameter in augmented Lagrangian approach","PCFieldSplitGKBNu",jac->gkbnu,&jac->gkbnu,NULL);
1645: PetscOptionsInt("-pc_fieldsplit_gkb_maxit","Maximum allowed number of iterations","PCFieldSplitGKBMaxit",jac->gkbmaxit,&jac->gkbmaxit,NULL);
1646: PetscOptionsBool("-pc_fieldsplit_gkb_monitor","Prints number of GKB iterations and error","PCFieldSplitGKB",jac->gkbmonitor,&jac->gkbmonitor,NULL);
1647: }
1648: PetscOptionsTail();
1649: return 0;
1650: }
1652: /*------------------------------------------------------------------------------------*/
1654: static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1655: {
1656: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1657: PC_FieldSplitLink ilink,next = jac->head;
1658: char prefix[128];
1659: PetscInt i;
1661: if (jac->splitdefined) {
1662: PetscInfo(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1663: return 0;
1664: }
1665: for (i=0; i<n; i++) {
1668: }
1669: PetscNew(&ilink);
1670: if (splitname) {
1671: PetscStrallocpy(splitname,&ilink->splitname);
1672: } else {
1673: PetscMalloc1(3,&ilink->splitname);
1674: PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);
1675: }
1676: 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 */
1677: PetscMalloc1(n,&ilink->fields);
1678: PetscArraycpy(ilink->fields,fields,n);
1679: PetscMalloc1(n,&ilink->fields_col);
1680: PetscArraycpy(ilink->fields_col,fields_col,n);
1682: ilink->nfields = n;
1683: ilink->next = NULL;
1684: KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1685: KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);
1686: PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1687: KSPSetType(ilink->ksp,KSPPREONLY);
1688: PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);
1690: PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);
1691: KSPSetOptionsPrefix(ilink->ksp,prefix);
1693: if (!next) {
1694: jac->head = ilink;
1695: ilink->previous = NULL;
1696: } else {
1697: while (next->next) {
1698: next = next->next;
1699: }
1700: next->next = ilink;
1701: ilink->previous = next;
1702: }
1703: jac->nsplits++;
1704: return 0;
1705: }
1707: static PetscErrorCode PCFieldSplitSchurGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp)
1708: {
1709: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1711: *subksp = NULL;
1712: if (n) *n = 0;
1713: if (jac->type == PC_COMPOSITE_SCHUR) {
1714: PetscInt nn;
1718: nn = jac->nsplits + (jac->kspupper != jac->head->ksp ? 1 : 0);
1719: PetscMalloc1(nn,subksp);
1720: (*subksp)[0] = jac->head->ksp;
1721: (*subksp)[1] = jac->kspschur;
1722: if (jac->kspupper != jac->head->ksp) (*subksp)[2] = jac->kspupper;
1723: if (n) *n = nn;
1724: }
1725: return 0;
1726: }
1728: static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp)
1729: {
1730: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1733: PetscMalloc1(jac->nsplits,subksp);
1734: MatSchurComplementGetKSP(jac->schur,*subksp);
1736: (*subksp)[1] = jac->kspschur;
1737: if (n) *n = jac->nsplits;
1738: return 0;
1739: }
1741: static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp)
1742: {
1743: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1744: PetscInt cnt = 0;
1745: PC_FieldSplitLink ilink = jac->head;
1747: PetscMalloc1(jac->nsplits,subksp);
1748: while (ilink) {
1749: (*subksp)[cnt++] = ilink->ksp;
1750: ilink = ilink->next;
1751: }
1753: if (n) *n = jac->nsplits;
1754: return 0;
1755: }
1757: /*@C
1758: PCFieldSplitRestrictIS - Restricts the fieldsplit ISs to be within a given IS.
1760: Input Parameters:
1761: + pc - the preconditioner context
1762: - is - the index set that defines the indices to which the fieldsplit is to be restricted
1764: Level: advanced
1766: @*/
1767: PetscErrorCode PCFieldSplitRestrictIS(PC pc,IS isy)
1768: {
1771: PetscTryMethod(pc,"PCFieldSplitRestrictIS_C",(PC,IS),(pc,isy));
1772: return 0;
1773: }
1775: static PetscErrorCode PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy)
1776: {
1777: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1778: PC_FieldSplitLink ilink = jac->head, next;
1779: PetscInt localsize,size,sizez,i;
1780: const PetscInt *ind, *indz;
1781: PetscInt *indc, *indcz;
1782: PetscBool flg;
1784: ISGetLocalSize(isy,&localsize);
1785: MPI_Scan(&localsize,&size,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)isy));
1786: size -= localsize;
1787: while (ilink) {
1788: IS isrl,isr;
1789: PC subpc;
1790: ISEmbed(ilink->is, isy, PETSC_TRUE, &isrl);
1791: ISGetLocalSize(isrl,&localsize);
1792: PetscMalloc1(localsize,&indc);
1793: ISGetIndices(isrl,&ind);
1794: PetscArraycpy(indc,ind,localsize);
1795: ISRestoreIndices(isrl,&ind);
1796: ISDestroy(&isrl);
1797: for (i=0; i<localsize; i++) *(indc+i) += size;
1798: ISCreateGeneral(PetscObjectComm((PetscObject)isy),localsize,indc,PETSC_OWN_POINTER,&isr);
1799: PetscObjectReference((PetscObject)isr);
1800: ISDestroy(&ilink->is);
1801: ilink->is = isr;
1802: PetscObjectReference((PetscObject)isr);
1803: ISDestroy(&ilink->is_col);
1804: ilink->is_col = isr;
1805: ISDestroy(&isr);
1806: KSPGetPC(ilink->ksp, &subpc);
1807: PetscObjectTypeCompare((PetscObject)subpc,PCFIELDSPLIT,&flg);
1808: if (flg) {
1809: IS iszl,isz;
1810: MPI_Comm comm;
1811: ISGetLocalSize(ilink->is,&localsize);
1812: comm = PetscObjectComm((PetscObject)ilink->is);
1813: ISEmbed(isy, ilink->is, PETSC_TRUE, &iszl);
1814: MPI_Scan(&localsize,&sizez,1,MPIU_INT,MPI_SUM,comm);
1815: sizez -= localsize;
1816: ISGetLocalSize(iszl,&localsize);
1817: PetscMalloc1(localsize,&indcz);
1818: ISGetIndices(iszl,&indz);
1819: PetscArraycpy(indcz,indz,localsize);
1820: ISRestoreIndices(iszl,&indz);
1821: ISDestroy(&iszl);
1822: for (i=0; i<localsize; i++) *(indcz+i) += sizez;
1823: ISCreateGeneral(comm,localsize,indcz,PETSC_OWN_POINTER,&isz);
1824: PCFieldSplitRestrictIS(subpc,isz);
1825: ISDestroy(&isz);
1826: }
1827: next = ilink->next;
1828: ilink = next;
1829: }
1830: jac->isrestrict = PETSC_TRUE;
1831: return 0;
1832: }
1834: static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is)
1835: {
1836: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1837: PC_FieldSplitLink ilink, next = jac->head;
1838: char prefix[128];
1840: if (jac->splitdefined) {
1841: PetscInfo(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1842: return 0;
1843: }
1844: PetscNew(&ilink);
1845: if (splitname) {
1846: PetscStrallocpy(splitname,&ilink->splitname);
1847: } else {
1848: PetscMalloc1(8,&ilink->splitname);
1849: PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);
1850: }
1851: 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 */
1852: PetscObjectReference((PetscObject)is);
1853: ISDestroy(&ilink->is);
1854: ilink->is = is;
1855: PetscObjectReference((PetscObject)is);
1856: ISDestroy(&ilink->is_col);
1857: ilink->is_col = is;
1858: ilink->next = NULL;
1859: KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1860: KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);
1861: PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1862: KSPSetType(ilink->ksp,KSPPREONLY);
1863: PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);
1865: PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);
1866: KSPSetOptionsPrefix(ilink->ksp,prefix);
1868: if (!next) {
1869: jac->head = ilink;
1870: ilink->previous = NULL;
1871: } else {
1872: while (next->next) {
1873: next = next->next;
1874: }
1875: next->next = ilink;
1876: ilink->previous = next;
1877: }
1878: jac->nsplits++;
1879: return 0;
1880: }
1882: /*@C
1883: PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner
1885: Logically Collective on PC
1887: Input Parameters:
1888: + pc - the preconditioner context
1889: . splitname - name of this split, if NULL the number of the split is used
1890: . n - the number of fields in this split
1891: - fields - the fields in this split
1893: Level: intermediate
1895: Notes:
1896: Use PCFieldSplitSetIS() to set a completely general set of indices as a field.
1898: The PCFieldSplitSetFields() is for defining fields as strided blocks. For example, if the block
1899: 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
1900: 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x....
1901: where the numbered entries indicate what is in the field.
1903: This function is called once per split (it creates a new split each time). Solve options
1904: for this split will be available under the prefix -fieldsplit_SPLITNAME_.
1906: Developer Note: This routine does not actually create the IS representing the split, that is delayed
1907: until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be
1908: available when this routine is called.
1910: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS()
1912: @*/
1913: PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1914: {
1919: PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt*,const PetscInt*),(pc,splitname,n,fields,fields_col));
1920: return 0;
1921: }
1923: /*@
1924: PCFieldSplitSetDiagUseAmat - set flag indicating whether to extract diagonal blocks from Amat (rather than Pmat)
1926: Logically Collective on PC
1928: Input Parameters:
1929: + pc - the preconditioner object
1930: - flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from
1932: Options Database:
1933: . -pc_fieldsplit_diag_use_amat - use the Amat to provide the diagonal blocks
1935: Level: intermediate
1937: .seealso: PCFieldSplitGetDiagUseAmat(), PCFieldSplitSetOffDiagUseAmat(), PCFIELDSPLIT
1939: @*/
1940: PetscErrorCode PCFieldSplitSetDiagUseAmat(PC pc,PetscBool flg)
1941: {
1942: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1943: PetscBool isfs;
1946: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1948: jac->diag_use_amat = flg;
1949: return 0;
1950: }
1952: /*@
1953: PCFieldSplitGetDiagUseAmat - get the flag indicating whether to extract diagonal blocks from Amat (rather than Pmat)
1955: Logically Collective on PC
1957: Input Parameters:
1958: . pc - the preconditioner object
1960: Output Parameters:
1961: . flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from
1963: Level: intermediate
1965: .seealso: PCFieldSplitSetDiagUseAmat(), PCFieldSplitGetOffDiagUseAmat(), PCFIELDSPLIT
1967: @*/
1968: PetscErrorCode PCFieldSplitGetDiagUseAmat(PC pc,PetscBool *flg)
1969: {
1970: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1971: PetscBool isfs;
1975: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1977: *flg = jac->diag_use_amat;
1978: return 0;
1979: }
1981: /*@
1982: PCFieldSplitSetOffDiagUseAmat - set flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat)
1984: Logically Collective on PC
1986: Input Parameters:
1987: + pc - the preconditioner object
1988: - flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from
1990: Options Database:
1991: . -pc_fieldsplit_off_diag_use_amat <bool> - use the Amat to extract the off-diagonal blocks
1993: Level: intermediate
1995: .seealso: PCFieldSplitGetOffDiagUseAmat(), PCFieldSplitSetDiagUseAmat(), PCFIELDSPLIT
1997: @*/
1998: PetscErrorCode PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg)
1999: {
2000: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2001: PetscBool isfs;
2004: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2006: jac->offdiag_use_amat = flg;
2007: return 0;
2008: }
2010: /*@
2011: PCFieldSplitGetOffDiagUseAmat - get the flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat)
2013: Logically Collective on PC
2015: Input Parameters:
2016: . pc - the preconditioner object
2018: Output Parameters:
2019: . 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;
2033: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2035: *flg = jac->offdiag_use_amat;
2036: return 0;
2037: }
2039: /*@C
2040: PCFieldSplitSetIS - Sets the exact elements for field
2042: Logically Collective on PC
2044: Input Parameters:
2045: + pc - the preconditioner context
2046: . splitname - name of this split, if NULL the number of the split is used
2047: - is - the index set that defines the vector elements in this field
2049: Notes:
2050: Use PCFieldSplitSetFields(), for fields defined by strided types.
2052: This function is called once per split (it creates a new split each time). Solve options
2053: for this split will be available under the prefix -fieldsplit_SPLITNAME_.
2055: Level: intermediate
2057: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize()
2059: @*/
2060: PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is)
2061: {
2065: PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));
2066: return 0;
2067: }
2069: /*@C
2070: PCFieldSplitGetIS - Retrieves the elements for a field as an IS
2072: Logically Collective on PC
2074: Input Parameters:
2075: + pc - the preconditioner context
2076: - splitname - name of this split
2078: Output Parameter:
2079: - is - the index set that defines the vector elements in this field, or NULL if the field is not found
2081: Level: intermediate
2083: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS()
2085: @*/
2086: PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is)
2087: {
2091: {
2092: PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;
2093: PC_FieldSplitLink ilink = jac->head;
2094: PetscBool found;
2096: *is = NULL;
2097: while (ilink) {
2098: PetscStrcmp(ilink->splitname, splitname, &found);
2099: if (found) {
2100: *is = ilink->is;
2101: break;
2102: }
2103: ilink = ilink->next;
2104: }
2105: }
2106: return 0;
2107: }
2109: /*@C
2110: PCFieldSplitGetISByIndex - Retrieves the elements for a given index field as an IS
2112: Logically Collective on PC
2114: Input Parameters:
2115: + pc - the preconditioner context
2116: - index - index of this split
2118: Output Parameter:
2119: - is - the index set that defines the vector elements in this field
2121: Level: intermediate
2123: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitGetIS(), PCFieldSplitSetIS()
2125: @*/
2126: PetscErrorCode PCFieldSplitGetISByIndex(PC pc,PetscInt index,IS *is)
2127: {
2131: {
2132: PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;
2133: PC_FieldSplitLink ilink = jac->head;
2134: PetscInt i = 0;
2137: while (i < index) {
2138: ilink = ilink->next;
2139: ++i;
2140: }
2141: PCFieldSplitGetIS(pc,ilink->splitname,is);
2142: }
2143: return 0;
2144: }
2146: /*@
2147: PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the
2148: fieldsplit preconditioner. If not set the matrix block size is used.
2150: Logically Collective on PC
2152: Input Parameters:
2153: + pc - the preconditioner context
2154: - bs - the block size
2156: Level: intermediate
2158: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields()
2160: @*/
2161: PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs)
2162: {
2165: PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));
2166: return 0;
2167: }
2169: /*@C
2170: PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits
2172: Collective on KSP
2174: Input Parameter:
2175: . pc - the preconditioner context
2177: Output Parameters:
2178: + n - the number of splits
2179: - subksp - the array of KSP contexts
2181: Note:
2182: After PCFieldSplitGetSubKSP() the array of KSPs is to be freed by the user with PetscFree()
2183: (not the KSP just the array that contains them).
2185: You must call PCSetUp() before calling PCFieldSplitGetSubKSP().
2187: If the fieldsplit is of type PC_COMPOSITE_SCHUR, it returns the KSP object used inside the
2188: Schur complement and the KSP object used to iterate over the Schur complement.
2189: To access all the KSP objects used in PC_COMPOSITE_SCHUR, use PCFieldSplitSchurGetSubKSP().
2191: If the fieldsplit is of type PC_COMPOSITE_GKB, it returns the KSP object used to solve the
2192: inner linear system defined by the matrix H in each loop.
2194: Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs.
2195: You can call PCFieldSplitGetSubKSP(pc,n,PETSC_NULL_KSP,ierr) to determine how large the
2196: KSP array must be.
2198: Level: advanced
2200: .seealso: PCFIELDSPLIT
2201: @*/
2202: PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
2203: {
2206: PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));
2207: return 0;
2208: }
2210: /*@C
2211: PCFieldSplitSchurGetSubKSP - Gets the KSP contexts used inside the Schur complement based PCFIELDSPLIT
2213: Collective on KSP
2215: Input Parameter:
2216: . pc - the preconditioner context
2218: Output Parameters:
2219: + n - the number of splits
2220: - subksp - the array of KSP contexts
2222: Note:
2223: After PCFieldSplitSchurGetSubKSP() the array of KSPs is to be freed by the user with PetscFree()
2224: (not the KSP just the array that contains them).
2226: You must call PCSetUp() before calling PCFieldSplitSchurGetSubKSP().
2228: If the fieldsplit type is of type PC_COMPOSITE_SCHUR, it returns (in order)
2229: - the KSP used for the (1,1) block
2230: - the KSP used for the Schur complement (not the one used for the interior Schur solver)
2231: - the KSP used for the (1,1) block in the upper triangular factor (if different from that of the (1,1) block).
2233: It returns a null array if the fieldsplit is not of type PC_COMPOSITE_SCHUR; in this case, you should use PCFieldSplitGetSubKSP().
2235: Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs.
2236: You can call PCFieldSplitSchurGetSubKSP(pc,n,PETSC_NULL_KSP,ierr) to determine how large the
2237: KSP array must be.
2239: Level: advanced
2241: .seealso: PCFIELDSPLIT
2242: @*/
2243: PetscErrorCode PCFieldSplitSchurGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
2244: {
2247: PetscUseMethod(pc,"PCFieldSplitSchurGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));
2248: return 0;
2249: }
2251: /*@
2252: PCFieldSplitSetSchurPre - Indicates from what operator the preconditioner is constructucted for the Schur complement.
2253: The default is the A11 matrix.
2255: Collective on PC
2257: Input Parameters:
2258: + pc - the preconditioner context
2259: . 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
2260: PC_FIELDSPLIT_SCHUR_PRE_SELFP, and PC_FIELDSPLIT_SCHUR_PRE_FULL
2261: - userpre - matrix to use for preconditioning, or NULL
2263: Options Database:
2264: + -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11. See notes for meaning of various arguments
2265: - -fieldsplit_1_pc_type <pctype> - the preconditioner algorithm that is used to construct the preconditioner from the operator
2267: Notes:
2268: $ If ptype is
2269: $ a11 - the preconditioner for the Schur complement is generated from the block diagonal part of the preconditioner
2270: $ matrix associated with the Schur complement (i.e. A11), not the Schur complement matrix
2271: $ self - the preconditioner for the Schur complement is generated from the symbolic representation of the Schur complement matrix:
2272: $ The only preconditioner that currently works with this symbolic respresentation matrix object is the PCLSC
2273: $ preconditioner
2274: $ user - the preconditioner for the Schur complement is generated from the user provided matrix (pre argument
2275: $ to this function).
2276: $ selfp - the preconditioning for the Schur complement is generated from an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01
2277: $ This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be
2278: $ lumped before extracting the diagonal using the additional option -fieldsplit_1_mat_schur_complement_ainv_type lump
2279: $ full - the preconditioner for the Schur complement is generated from the exact Schur complement matrix representation computed internally by PCFIELDSPLIT (this is expensive)
2280: $ useful mostly as a test that the Schur complement approach can work for your problem
2282: When solving a saddle point problem, where the A11 block is identically zero, using a11 as the ptype only makes sense
2283: with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and
2284: -fieldsplit_1_pc_type lsc which uses the least squares commutator to compute a preconditioner for the Schur complement.
2286: Level: intermediate
2288: .seealso: PCFieldSplitGetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType,
2289: MatSchurComplementSetAinvType(), PCLSC
2291: @*/
2292: PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
2293: {
2295: PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));
2296: return 0;
2297: }
2299: PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) {return PCFieldSplitSetSchurPre(pc,ptype,pre);} /* Deprecated name */
2301: /*@
2302: PCFieldSplitGetSchurPre - For Schur complement fieldsplit, determine how the Schur complement will be
2303: preconditioned. See PCFieldSplitSetSchurPre() for details.
2305: Logically Collective on PC
2307: Input Parameter:
2308: . pc - the preconditioner context
2310: Output Parameters:
2311: + ptype - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11, PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_PRE_USER
2312: - userpre - matrix to use for preconditioning (with PC_FIELDSPLIT_PRE_USER), or NULL
2314: Level: intermediate
2316: .seealso: PCFieldSplitSetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC
2318: @*/
2319: PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
2320: {
2322: PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre));
2323: return 0;
2324: }
2326: /*@
2327: PCFieldSplitSchurGetS - extract the MatSchurComplement object used by this PC in case it needs to be configured separately
2329: Not collective
2331: Input Parameter:
2332: . pc - the preconditioner context
2334: Output Parameter:
2335: . S - the Schur complement matrix
2337: Notes:
2338: This matrix should not be destroyed using MatDestroy(); rather, use PCFieldSplitSchurRestoreS().
2340: Level: advanced
2342: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurRestoreS()
2344: @*/
2345: PetscErrorCode PCFieldSplitSchurGetS(PC pc,Mat *S)
2346: {
2347: const char* t;
2348: PetscBool isfs;
2349: PC_FieldSplit *jac;
2352: PetscObjectGetType((PetscObject)pc,&t);
2353: PetscStrcmp(t,PCFIELDSPLIT,&isfs);
2355: jac = (PC_FieldSplit*)pc->data;
2357: if (S) *S = jac->schur;
2358: return 0;
2359: }
2361: /*@
2362: PCFieldSplitSchurRestoreS - restores the MatSchurComplement object used by this PC
2364: Not collective
2366: Input Parameters:
2367: + pc - the preconditioner context
2368: - S - the Schur complement matrix
2370: Level: advanced
2372: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurGetS()
2374: @*/
2375: PetscErrorCode PCFieldSplitSchurRestoreS(PC pc,Mat *S)
2376: {
2377: const char* t;
2378: PetscBool isfs;
2379: PC_FieldSplit *jac;
2382: PetscObjectGetType((PetscObject)pc,&t);
2383: PetscStrcmp(t,PCFIELDSPLIT,&isfs);
2385: jac = (PC_FieldSplit*)pc->data;
2388: return 0;
2389: }
2391: static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
2392: {
2393: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2395: jac->schurpre = ptype;
2396: if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) {
2397: MatDestroy(&jac->schur_user);
2398: jac->schur_user = pre;
2399: PetscObjectReference((PetscObject)jac->schur_user);
2400: }
2401: return 0;
2402: }
2404: static PetscErrorCode PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
2405: {
2406: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2408: *ptype = jac->schurpre;
2409: *pre = jac->schur_user;
2410: return 0;
2411: }
2413: /*@
2414: PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain in the preconditioner
2416: Collective on PC
2418: Input Parameters:
2419: + pc - the preconditioner context
2420: - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default
2422: Options Database:
2423: . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> - default is full
2425: Level: intermediate
2427: Notes:
2428: The FULL factorization is
2430: .vb
2431: (A B) = (1 0) (A 0) (1 Ainv*B) = L D U
2432: (C E) (C*Ainv 1) (0 S) (0 1)
2433: .vb
2434: 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,
2435: 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().
2437: If A and S are solved exactly
2438: .vb
2439: *) FULL factorization is a direct solver.
2440: *) 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.
2441: *) 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.
2442: .ve
2444: If the iteration count is very low, consider using KSPFGMRES or KSPGCR which can use one less preconditioner
2445: application in this case. Note that the preconditioned operator may be highly non-normal, so such fast convergence may not be observed in practice.
2447: For symmetric problems in which A is positive definite and S is negative definite, DIAG can be used with KSPMINRES.
2449: 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).
2451: References:
2452: + * - Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000).
2453: - * - Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001).
2455: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCFieldSplitSetSchurScale()
2456: @*/
2457: PetscErrorCode PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype)
2458: {
2460: PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));
2461: return 0;
2462: }
2464: static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype)
2465: {
2466: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2468: jac->schurfactorization = ftype;
2469: return 0;
2470: }
2472: /*@
2473: PCFieldSplitSetSchurScale - Controls the sign flip of S for PC_FIELDSPLIT_SCHUR_FACT_DIAG.
2475: Collective on PC
2477: Input Parameters:
2478: + pc - the preconditioner context
2479: - scale - scaling factor for the Schur complement
2481: Options Database:
2482: . -pc_fieldsplit_schur_scale - default is -1.0
2484: Level: intermediate
2486: .seealso: PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurFactType, PCFieldSplitSetSchurScale()
2487: @*/
2488: PetscErrorCode PCFieldSplitSetSchurScale(PC pc,PetscScalar scale)
2489: {
2492: PetscTryMethod(pc,"PCFieldSplitSetSchurScale_C",(PC,PetscScalar),(pc,scale));
2493: return 0;
2494: }
2496: static PetscErrorCode PCFieldSplitSetSchurScale_FieldSplit(PC pc,PetscScalar scale)
2497: {
2498: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2500: jac->schurscale = scale;
2501: return 0;
2502: }
2504: /*@C
2505: PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement
2507: Collective on KSP
2509: Input Parameter:
2510: . pc - the preconditioner context
2512: Output Parameters:
2513: + A00 - the (0,0) block
2514: . A01 - the (0,1) block
2515: . A10 - the (1,0) block
2516: - A11 - the (1,1) block
2518: Level: advanced
2520: .seealso: PCFIELDSPLIT
2521: @*/
2522: PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11)
2523: {
2524: PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;
2528: if (A00) *A00 = jac->pmat[0];
2529: if (A01) *A01 = jac->B;
2530: if (A10) *A10 = jac->C;
2531: if (A11) *A11 = jac->pmat[1];
2532: return 0;
2533: }
2535: /*@
2536: PCFieldSplitSetGKBTol - Sets the solver tolerance for the generalized Golub-Kahan bidiagonalization preconditioner.
2538: Collective on PC
2540: Notes:
2541: The generalized GKB algorithm uses a lower bound estimate of the error in energy norm as stopping criterion.
2542: It stops once the lower bound estimate undershoots the required solver tolerance. Although the actual error might be bigger than
2543: this estimate, the stopping criterion is satisfactory in practical cases [A13].
2545: [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013.
2547: Input Parameters:
2548: + pc - the preconditioner context
2549: - tolerance - the solver tolerance
2551: Options Database:
2552: . -pc_fieldsplit_gkb_tol - default is 1e-5
2554: Level: intermediate
2556: .seealso: PCFIELDSPLIT, PCFieldSplitSetGKBDelay(), PCFieldSplitSetGKBNu(), PCFieldSplitSetGKBMaxit()
2557: @*/
2558: PetscErrorCode PCFieldSplitSetGKBTol(PC pc,PetscReal tolerance)
2559: {
2562: PetscTryMethod(pc,"PCFieldSplitSetGKBTol_C",(PC,PetscReal),(pc,tolerance));
2563: return 0;
2564: }
2566: static PetscErrorCode PCFieldSplitSetGKBTol_FieldSplit(PC pc,PetscReal tolerance)
2567: {
2568: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2570: jac->gkbtol = tolerance;
2571: return 0;
2572: }
2574: /*@
2575: PCFieldSplitSetGKBMaxit - Sets the maximum number of iterations for the generalized Golub-Kahan bidiagonalization
2576: preconditioner.
2578: Collective on PC
2580: Input Parameters:
2581: + pc - the preconditioner context
2582: - maxit - the maximum number of iterations
2584: Options Database:
2585: . -pc_fieldsplit_gkb_maxit - default is 100
2587: Level: intermediate
2589: .seealso: PCFIELDSPLIT, PCFieldSplitSetGKBDelay(), PCFieldSplitSetGKBTol(), PCFieldSplitSetGKBNu()
2590: @*/
2591: PetscErrorCode PCFieldSplitSetGKBMaxit(PC pc,PetscInt maxit)
2592: {
2595: PetscTryMethod(pc,"PCFieldSplitSetGKBMaxit_C",(PC,PetscInt),(pc,maxit));
2596: return 0;
2597: }
2599: static PetscErrorCode PCFieldSplitSetGKBMaxit_FieldSplit(PC pc,PetscInt maxit)
2600: {
2601: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2603: jac->gkbmaxit = maxit;
2604: return 0;
2605: }
2607: /*@
2608: PCFieldSplitSetGKBDelay - Sets the delay in the lower bound error estimate in the generalized Golub-Kahan bidiagonalization
2609: preconditioner.
2611: Collective on PC
2613: Notes:
2614: 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
2615: is expressed as a truncated sum. The error at iteration k can only be measured at iteration (k + delay), and thus the algorithm needs
2616: 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
2618: [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013.
2620: Input Parameters:
2621: + pc - the preconditioner context
2622: - delay - the delay window in the lower bound estimate
2624: Options Database:
2625: . -pc_fieldsplit_gkb_delay - default is 5
2627: Level: intermediate
2629: .seealso: PCFIELDSPLIT, PCFieldSplitSetGKBNu(), PCFieldSplitSetGKBTol(), PCFieldSplitSetGKBMaxit()
2630: @*/
2631: PetscErrorCode PCFieldSplitSetGKBDelay(PC pc,PetscInt delay)
2632: {
2635: PetscTryMethod(pc,"PCFieldSplitSetGKBDelay_C",(PC,PetscInt),(pc,delay));
2636: return 0;
2637: }
2639: static PetscErrorCode PCFieldSplitSetGKBDelay_FieldSplit(PC pc,PetscInt delay)
2640: {
2641: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2643: jac->gkbdelay = delay;
2644: return 0;
2645: }
2647: /*@
2648: 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.
2650: Collective on PC
2652: Notes:
2653: 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,
2654: 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
2655: necessary to find a good balance in between the convergence of the inner and outer loop.
2657: 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.
2659: [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013.
2661: Input Parameters:
2662: + pc - the preconditioner context
2663: - nu - the shift parameter
2665: Options Database:
2666: . -pc_fieldsplit_gkb_nu - default is 1
2668: Level: intermediate
2670: .seealso: PCFIELDSPLIT, PCFieldSplitSetGKBDelay(), PCFieldSplitSetGKBTol(), PCFieldSplitSetGKBMaxit()
2671: @*/
2672: PetscErrorCode PCFieldSplitSetGKBNu(PC pc,PetscReal nu)
2673: {
2676: PetscTryMethod(pc,"PCFieldSplitSetGKBNu_C",(PC,PetscReal),(pc,nu));
2677: return 0;
2678: }
2680: static PetscErrorCode PCFieldSplitSetGKBNu_FieldSplit(PC pc,PetscReal nu)
2681: {
2682: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2684: jac->gkbnu = nu;
2685: return 0;
2686: }
2688: static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type)
2689: {
2690: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2692: jac->type = type;
2694: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",0);
2695: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",0);
2696: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",0);
2697: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);
2698: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",0);
2699: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBTol_C",0);
2700: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBMaxit_C",0);
2701: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBNu_C",0);
2702: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBDelay_C",0);
2704: if (type == PC_COMPOSITE_SCHUR) {
2705: pc->ops->apply = PCApply_FieldSplit_Schur;
2706: pc->ops->view = PCView_FieldSplit_Schur;
2708: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);
2709: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit);
2710: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit);
2711: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);
2712: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",PCFieldSplitSetSchurScale_FieldSplit);
2713: } else if (type == PC_COMPOSITE_GKB) {
2714: pc->ops->apply = PCApply_FieldSplit_GKB;
2715: pc->ops->view = PCView_FieldSplit_GKB;
2717: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2718: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBTol_C",PCFieldSplitSetGKBTol_FieldSplit);
2719: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBMaxit_C",PCFieldSplitSetGKBMaxit_FieldSplit);
2720: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBNu_C",PCFieldSplitSetGKBNu_FieldSplit);
2721: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBDelay_C",PCFieldSplitSetGKBDelay_FieldSplit);
2722: } else {
2723: pc->ops->apply = PCApply_FieldSplit;
2724: pc->ops->view = PCView_FieldSplit;
2726: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2727: }
2728: return 0;
2729: }
2731: static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs)
2732: {
2733: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2737: jac->bs = bs;
2738: return 0;
2739: }
2741: static PetscErrorCode PCSetCoordinates_FieldSplit(PC pc, PetscInt dim, PetscInt nloc, PetscReal coords[])
2742: {
2743: PC_FieldSplit * jac = (PC_FieldSplit*)pc->data;
2744: PC_FieldSplitLink ilink_current = jac->head;
2745: PetscInt ii;
2746: IS is_owned;
2748: jac->coordinates_set = PETSC_TRUE; // Internal flag
2749: MatGetOwnershipIS(pc->mat,&is_owned,PETSC_NULL);
2751: ii=0;
2752: while (ilink_current) {
2753: // For each IS, embed it to get local coords indces
2754: IS is_coords;
2755: PetscInt ndofs_block;
2756: const PetscInt *block_dofs_enumeration; // Numbering of the dofs relevant to the current block
2758: // Setting drop to true for safety. It should make no difference.
2759: ISEmbed(ilink_current->is, is_owned, PETSC_TRUE, &is_coords);
2760: ISGetLocalSize(is_coords, &ndofs_block);
2761: ISGetIndices(is_coords, &block_dofs_enumeration);
2763: // Allocate coordinates vector and set it directly
2764: PetscMalloc1(ndofs_block * dim, &(ilink_current->coords));
2765: for (PetscInt dof=0;dof<ndofs_block;++dof) {
2766: for (PetscInt d=0;d<dim;++d) {
2767: (ilink_current->coords)[dim*dof + d] = coords[dim * block_dofs_enumeration[dof] + d];
2768: }
2769: }
2770: ilink_current->dim = dim;
2771: ilink_current->ndofs = ndofs_block;
2772: ISRestoreIndices(is_coords, &block_dofs_enumeration);
2773: ISDestroy(&is_coords);
2774: ilink_current = ilink_current->next;
2775: ++ii;
2776: }
2777: ISDestroy(&is_owned);
2778: return 0;
2779: }
2781: /*@
2782: PCFieldSplitSetType - Sets the type of fieldsplit preconditioner.
2784: Collective on PC
2786: Input Parameters:
2787: + pc - the preconditioner context
2788: - type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR
2790: Options Database Key:
2791: . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type
2793: Level: Intermediate
2795: .seealso: PCCompositeSetType()
2797: @*/
2798: PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type)
2799: {
2801: PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));
2802: return 0;
2803: }
2805: /*@
2806: PCFieldSplitGetType - Gets the type of fieldsplit preconditioner.
2808: Not collective
2810: Input Parameter:
2811: . pc - the preconditioner context
2813: Output Parameter:
2814: . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR
2816: Level: Intermediate
2818: .seealso: PCCompositeSetType()
2819: @*/
2820: PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type)
2821: {
2822: PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;
2826: *type = jac->type;
2827: return 0;
2828: }
2830: /*@
2831: PCFieldSplitSetDMSplits - Flags whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible.
2833: Logically Collective
2835: Input Parameters:
2836: + pc - the preconditioner context
2837: - flg - boolean indicating whether to use field splits defined by the DM
2839: Options Database Key:
2840: . -pc_fieldsplit_dm_splits <bool> - use the field splits defined by the DM
2842: Level: Intermediate
2844: .seealso: PCFieldSplitGetDMSplits()
2846: @*/
2847: PetscErrorCode PCFieldSplitSetDMSplits(PC pc,PetscBool flg)
2848: {
2849: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2850: PetscBool isfs;
2854: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2855: if (isfs) {
2856: jac->dm_splits = flg;
2857: }
2858: return 0;
2859: }
2861: /*@
2862: PCFieldSplitGetDMSplits - Returns flag indicating whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible.
2864: Logically Collective
2866: Input Parameter:
2867: . pc - the preconditioner context
2869: Output Parameter:
2870: . flg - boolean indicating whether to use field splits defined by the DM
2872: Level: Intermediate
2874: .seealso: PCFieldSplitSetDMSplits()
2876: @*/
2877: PetscErrorCode PCFieldSplitGetDMSplits(PC pc,PetscBool* flg)
2878: {
2879: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2880: PetscBool isfs;
2884: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2885: if (isfs) {
2886: if (flg) *flg = jac->dm_splits;
2887: }
2888: return 0;
2889: }
2891: /*@
2892: PCFieldSplitGetDetectSaddlePoint - Returns flag indicating whether PCFieldSplit will attempt to automatically determine fields based on zero diagonal entries.
2894: Logically Collective
2896: Input Parameter:
2897: . pc - the preconditioner context
2899: Output Parameter:
2900: . flg - boolean indicating whether to detect fields or not
2902: Level: Intermediate
2904: .seealso: PCFIELDSPLIT, PCFieldSplitSetDetectSaddlePoint()
2906: @*/
2907: PetscErrorCode PCFieldSplitGetDetectSaddlePoint(PC pc,PetscBool *flg)
2908: {
2909: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2911: *flg = jac->detect;
2912: return 0;
2913: }
2915: /*@
2916: PCFieldSplitSetDetectSaddlePoint - Sets flag indicating whether PCFieldSplit will attempt to automatically determine fields based on zero diagonal entries.
2918: Logically Collective
2920: Input Parameter:
2921: . pc - the preconditioner context
2923: Output Parameter:
2924: . flg - boolean indicating whether to detect fields or not
2926: Options Database Key:
2927: . -pc_fieldsplit_detect_saddle_point <bool> - detect and use the saddle point
2929: Notes:
2930: Also sets the split type to PC_COMPOSITE_SCHUR (see PCFieldSplitSetType()) and the Schur preconditioner type to PC_FIELDSPLIT_SCHUR_PRE_SELF (see PCFieldSplitSetSchurPre()).
2932: Level: Intermediate
2934: .seealso: PCFIELDSPLIT, PCFieldSplitSetDetectSaddlePoint(), PCFieldSplitSetType(), PCFieldSplitSetSchurPre()
2936: @*/
2937: PetscErrorCode PCFieldSplitSetDetectSaddlePoint(PC pc,PetscBool flg)
2938: {
2939: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2941: jac->detect = flg;
2942: if (jac->detect) {
2943: PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);
2944: PCFieldSplitSetSchurPre(pc,PC_FIELDSPLIT_SCHUR_PRE_SELF,NULL);
2945: }
2946: return 0;
2947: }
2949: /* -------------------------------------------------------------------------------------*/
2950: /*MC
2951: PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual
2952: fields or groups of fields. See the users manual section "Solving Block Matrices" for more details.
2954: To set options on the solvers for each block append -fieldsplit_ to all the PC
2955: options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1
2957: To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP()
2958: and set the options directly on the resulting KSP object
2960: Level: intermediate
2962: Options Database Keys:
2963: + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split
2964: . -pc_fieldsplit_default - automatically add any fields to additional splits that have not
2965: been supplied explicitly by -pc_fieldsplit_%d_fields
2966: . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields)
2967: . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur,gkb> - type of relaxation or factorization splitting
2968: . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11; see PCFieldSplitSetSchurPre()
2969: . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> - set factorization type when using -pc_fieldsplit_type schur; see PCFieldSplitSetSchurFactType()
2970: - -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero diagonal and uses Schur complement with no preconditioner as the solver
2972: Options prefixes for inner solvers when using the Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ .
2973: For all other solvers they are -fieldsplit_%d_ for the dth field; use -fieldsplit_ for all fields.
2974: The options prefix for the inner solver when using the Golub-Kahan biadiagonalization preconditioner is -fieldsplit_0_
2976: Notes:
2977: Use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS()
2978: to define a field by an arbitrary collection of entries.
2980: If no fields are set the default is used. The fields are defined by entries strided by bs,
2981: beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(),
2982: if this is not called the block size defaults to the blocksize of the second matrix passed
2983: to KSPSetOperators()/PCSetOperators().
2985: $ For the Schur complement preconditioner if J = [ A00 A01]
2986: $ [ A10 A11]
2987: $ the preconditioner using full factorization is
2988: $ [ I -ksp(A00) A01] [ inv(A00) 0 ] [ I 0]
2989: $ [ 0 I ] [ 0 ksp(S)] [ -A10 ksp(A00) I]
2990: where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. S is the Schur complement
2991: $ S = A11 - A10 ksp(A00) A01
2992: 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
2993: in providing the SECOND split or 1 if not give). For PCFieldSplitGetSubKSP() when field number is 0,
2994: it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default
2995: A11 is used to construct a preconditioner for S, use PCFieldSplitSetSchurPre() for all the possible ways to construct the preconditioner for S.
2997: The factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above,
2998: diag gives
2999: $ [ inv(A00) 0 ]
3000: $ [ 0 -ksp(S)]
3001: 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
3002: can be turned off with PCFieldSplitSetSchurScale() or by command line -pc_fieldsplit_schur_scale 1.0. The lower factorization is the inverse of
3003: $ [ A00 0]
3004: $ [ A10 S]
3005: where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of
3006: $ [ A00 A01]
3007: $ [ 0 S ]
3008: where again the inverses of A00 and S are applied using KSPs.
3010: If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS
3011: is used automatically for a second block.
3013: The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1.
3014: Generally it should be used with the AIJ format.
3016: The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see,
3017: for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT
3018: inside a smoother resulting in "Distributive Smoothers".
3020: References:
3021: . * - A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier-Stokes equations
3022: Howard Elman, V.E. Howle, John Shadid, Robert Shuttleworth, Ray Tuminaro, Journal of Computational Physics 227 (2008) 1790--1808
3023: http://chess.cs.umd.edu/~elman/papers/tax.pdf
3025: The Constrained Pressure Preconditioner (CPR) can be implemented using PCCOMPOSITE with PCGALERKIN. CPR first solves an R A P subsystem, updates the
3026: residual on all variables (PCCompositeSetType(pc,PC_COMPOSITE_MULTIPLICATIVE)), and then applies a simple ILU like preconditioner on all the variables.
3028: The generalized Golub-Kahan bidiagonalization preconditioner (gkb) can be applied to symmetric 2x2 block matrices of the shape
3029: $ [ A00 A01]
3030: $ [ A01' 0 ]
3031: 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'.
3032: 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_.
3034: . * - [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013.
3036: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCLSC,
3037: PCFieldSplitGetSubKSP(), PCFieldSplitSchurGetSubKSP(), PCFieldSplitSetFields(),
3038: PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre(), PCFieldSplitSetSchurFactType(),
3039: MatSchurComplementSetAinvType(), PCFieldSplitSetSchurScale(), PCFieldSplitSetDetectSaddlePoint()
3040: M*/
3042: PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc)
3043: {
3044: PC_FieldSplit *jac;
3046: PetscNewLog(pc,&jac);
3048: jac->bs = -1;
3049: jac->nsplits = 0;
3050: jac->type = PC_COMPOSITE_MULTIPLICATIVE;
3051: jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */
3052: jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL;
3053: jac->schurscale = -1.0;
3054: jac->dm_splits = PETSC_TRUE;
3055: jac->detect = PETSC_FALSE;
3056: jac->gkbtol = 1e-5;
3057: jac->gkbdelay = 5;
3058: jac->gkbnu = 1;
3059: jac->gkbmaxit = 100;
3060: jac->gkbmonitor = PETSC_FALSE;
3061: jac->coordinates_set = PETSC_FALSE;
3063: pc->data = (void*)jac;
3065: pc->ops->apply = PCApply_FieldSplit;
3066: pc->ops->applytranspose = PCApplyTranspose_FieldSplit;
3067: pc->ops->setup = PCSetUp_FieldSplit;
3068: pc->ops->reset = PCReset_FieldSplit;
3069: pc->ops->destroy = PCDestroy_FieldSplit;
3070: pc->ops->setfromoptions = PCSetFromOptions_FieldSplit;
3071: pc->ops->view = PCView_FieldSplit;
3072: pc->ops->applyrichardson = NULL;
3074: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurGetSubKSP_C",PCFieldSplitSchurGetSubKSP_FieldSplit);
3075: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
3076: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);
3077: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);
3078: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);
3079: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);
3080: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",PCFieldSplitRestrictIS_FieldSplit);
3081: PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_FieldSplit);
3082: return 0;
3083: }