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;
21: };
23: typedef struct {
24: PCCompositeType type;
25: PetscBool defaultsplit; /* Flag for a system with a set of 'k' scalar fields with the same layout (and bs = k) */
26: PetscBool splitdefined; /* Flag is set after the splits have been defined, to prevent more splits from being added */
27: PetscInt bs; /* Block size for IS and Mat structures */
28: PetscInt nsplits; /* Number of field divisions defined */
29: Vec *x,*y,w1,w2;
30: Mat *mat; /* The diagonal block for each split */
31: Mat *pmat; /* The preconditioning diagonal block for each split */
32: Mat *Afield; /* The rows of the matrix associated with each split */
33: PetscBool issetup;
35: /* Only used when Schur complement preconditioning is used */
36: Mat B; /* The (0,1) block */
37: Mat C; /* The (1,0) block */
38: Mat schur; /* The Schur complement S = A11 - A10 A00^{-1} A01, the KSP here, kspinner, is H_1 in [El08] */
39: Mat schurp; /* Assembled approximation to S built by MatSchurComplement to be used as a preconditioning matrix when solving with S */
40: Mat schur_user; /* User-provided preconditioning matrix for the Schur complement */
41: PCFieldSplitSchurPreType schurpre; /* Determines which preconditioning matrix is used for the Schur complement */
42: PCFieldSplitSchurFactType schurfactorization;
43: KSP kspschur; /* The solver for S */
44: KSP kspupper; /* The solver for A in the upper diagonal part of the factorization (H_2 in [El08]) */
45: PetscScalar schurscale; /* Scaling factor for the Schur complement solution with DIAG factorization */
47: /* Only used when Golub-Kahan bidiagonalization preconditioning is used */
48: Mat H; /* The modified matrix H = A00 + nu*A01*A01' */
49: PetscReal gkbtol; /* Stopping tolerance for lower bound estimate */
50: PetscInt gkbdelay; /* The delay window for the stopping criterion */
51: PetscReal gkbnu; /* Parameter for augmented Lagrangian H = A + nu*A01*A01' */
52: PetscInt gkbmaxit; /* Maximum number of iterations for outer loop */
53: PetscBool gkbmonitor; /* Monitor for gkb iterations and the lower bound error */
54: PetscViewer gkbviewer; /* Viewer context for gkbmonitor */
55: Vec u,v,d,Hu; /* Work vectors for the GKB algorithm */
56: PetscScalar *vecz; /* Contains intermediate values, eg for lower bound */
58: PC_FieldSplitLink head;
59: PetscBool isrestrict; /* indicates PCFieldSplitRestrictIS() has been last called on this object, hack */
60: PetscBool suboptionsset; /* Indicates that the KSPSetFromOptions() has been called on the sub-KSPs */
61: PetscBool dm_splits; /* Whether to use DMCreateFieldDecomposition() whenever possible */
62: PetscBool diag_use_amat; /* Whether to extract diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */
63: PetscBool offdiag_use_amat; /* Whether to extract off-diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */
64: PetscBool detect; /* Whether to form 2-way split by finding zero diagonal entries */
65: } PC_FieldSplit;
67: /*
68: Notes:
69: there is no particular reason that pmat, x, and y are stored as arrays in PC_FieldSplit instead of
70: inside PC_FieldSplitLink, just historical. If you want to be able to add new fields after already using the
71: PC you could change this.
72: */
74: /* This helper is so that setting a user-provided preconditioning matrix is orthogonal to choosing to use it. This way the
75: * application-provided FormJacobian can provide this matrix without interfering with the user's (command-line) choices. */
76: static Mat FieldSplitSchurPre(PC_FieldSplit *jac)
77: {
78: switch (jac->schurpre) {
79: case PC_FIELDSPLIT_SCHUR_PRE_SELF: return jac->schur;
80: case PC_FIELDSPLIT_SCHUR_PRE_SELFP: return jac->schurp;
81: case PC_FIELDSPLIT_SCHUR_PRE_A11: return jac->pmat[1];
82: case PC_FIELDSPLIT_SCHUR_PRE_FULL: /* We calculate this and store it in schur_user */
83: case PC_FIELDSPLIT_SCHUR_PRE_USER: /* Use a user-provided matrix if it is given, otherwise diagonal block */
84: default:
85: return jac->schur_user ? jac->schur_user : jac->pmat[1];
86: }
87: }
89: #include <petscdraw.h>
90: static PetscErrorCode PCView_FieldSplit(PC pc,PetscViewer viewer)
91: {
92: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
93: PetscErrorCode ierr;
94: PetscBool iascii,isdraw;
95: PetscInt i,j;
96: PC_FieldSplitLink ilink = jac->head;
99: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
100: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
101: if (iascii) {
102: if (jac->bs > 0) {
103: PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs);
104: } else {
105: PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D\n",PCCompositeTypes[jac->type],jac->nsplits);
106: }
107: if (pc->useAmat) {
108: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");
109: }
110: if (jac->diag_use_amat) {
111: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for diagonal blocks\n");
112: }
113: if (jac->offdiag_use_amat) {
114: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for off-diagonal blocks\n");
115: }
116: PetscViewerASCIIPrintf(viewer," Solver info for each split is in the following KSP objects:\n");
117: for (i=0; i<jac->nsplits; i++) {
118: if (ilink->fields) {
119: PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);
120: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
121: for (j=0; j<ilink->nfields; j++) {
122: if (j > 0) {
123: PetscViewerASCIIPrintf(viewer,",");
124: }
125: PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);
126: }
127: PetscViewerASCIIPrintf(viewer,"\n");
128: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
129: } else {
130: PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);
131: }
132: KSPView(ilink->ksp,viewer);
133: ilink = ilink->next;
134: }
135: }
137: if (isdraw) {
138: PetscDraw draw;
139: PetscReal x,y,w,wd;
141: PetscViewerDrawGetDraw(viewer,0,&draw);
142: PetscDrawGetCurrentPoint(draw,&x,&y);
143: w = 2*PetscMin(1.0 - x,x);
144: wd = w/(jac->nsplits + 1);
145: x = x - wd*(jac->nsplits-1)/2.0;
146: for (i=0; i<jac->nsplits; i++) {
147: PetscDrawPushCurrentPoint(draw,x,y);
148: KSPView(ilink->ksp,viewer);
149: PetscDrawPopCurrentPoint(draw);
150: x += wd;
151: ilink = ilink->next;
152: }
153: }
154: return(0);
155: }
157: static PetscErrorCode PCView_FieldSplit_Schur(PC pc,PetscViewer viewer)
158: {
159: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
160: PetscErrorCode ierr;
161: PetscBool iascii,isdraw;
162: PetscInt i,j;
163: PC_FieldSplitLink ilink = jac->head;
164: MatSchurComplementAinvType atype;
167: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
168: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
169: if (iascii) {
170: if (jac->bs > 0) {
171: PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, blocksize = %D, factorization %s\n",jac->bs,PCFieldSplitSchurFactTypes[jac->schurfactorization]);
172: } else {
173: PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, factorization %s\n",PCFieldSplitSchurFactTypes[jac->schurfactorization]);
174: }
175: if (pc->useAmat) {
176: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");
177: }
178: switch (jac->schurpre) {
179: case PC_FIELDSPLIT_SCHUR_PRE_SELF:
180: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from S itself\n");
181: break;
182: case PC_FIELDSPLIT_SCHUR_PRE_SELFP:
183: MatSchurComplementGetAinvType(jac->schur,&atype);
184: 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;
185: case PC_FIELDSPLIT_SCHUR_PRE_A11:
186: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n");
187: break;
188: case PC_FIELDSPLIT_SCHUR_PRE_FULL:
189: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from the exact Schur complement\n");
190: break;
191: case PC_FIELDSPLIT_SCHUR_PRE_USER:
192: if (jac->schur_user) {
193: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from user provided matrix\n");
194: } else {
195: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n");
196: }
197: break;
198: default:
199: SETERRQ1(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Invalid Schur preconditioning type: %d", jac->schurpre);
200: }
201: PetscViewerASCIIPrintf(viewer," Split info:\n");
202: PetscViewerASCIIPushTab(viewer);
203: for (i=0; i<jac->nsplits; i++) {
204: if (ilink->fields) {
205: PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);
206: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
207: for (j=0; j<ilink->nfields; j++) {
208: if (j > 0) {
209: PetscViewerASCIIPrintf(viewer,",");
210: }
211: PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);
212: }
213: PetscViewerASCIIPrintf(viewer,"\n");
214: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
215: } else {
216: PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);
217: }
218: ilink = ilink->next;
219: }
220: PetscViewerASCIIPrintf(viewer,"KSP solver for A00 block\n");
221: PetscViewerASCIIPushTab(viewer);
222: if (jac->head) {
223: KSPView(jac->head->ksp,viewer);
224: } else {PetscViewerASCIIPrintf(viewer," not yet available\n");}
225: PetscViewerASCIIPopTab(viewer);
226: if (jac->head && jac->kspupper != jac->head->ksp) {
227: PetscViewerASCIIPrintf(viewer,"KSP solver for upper A00 in upper triangular factor \n");
228: PetscViewerASCIIPushTab(viewer);
229: if (jac->kspupper) {KSPView(jac->kspupper,viewer);}
230: else {PetscViewerASCIIPrintf(viewer," not yet available\n");}
231: PetscViewerASCIIPopTab(viewer);
232: }
233: PetscViewerASCIIPrintf(viewer,"KSP solver for S = A11 - A10 inv(A00) A01 \n");
234: PetscViewerASCIIPushTab(viewer);
235: if (jac->kspschur) {
236: KSPView(jac->kspschur,viewer);
237: } else {
238: PetscViewerASCIIPrintf(viewer," not yet available\n");
239: }
240: PetscViewerASCIIPopTab(viewer);
241: PetscViewerASCIIPopTab(viewer);
242: } else if (isdraw && jac->head) {
243: PetscDraw draw;
244: PetscReal x,y,w,wd,h;
245: PetscInt cnt = 2;
246: char str[32];
248: PetscViewerDrawGetDraw(viewer,0,&draw);
249: PetscDrawGetCurrentPoint(draw,&x,&y);
250: if (jac->kspupper != jac->head->ksp) cnt++;
251: w = 2*PetscMin(1.0 - x,x);
252: wd = w/(cnt + 1);
254: PetscSNPrintf(str,32,"Schur fact. %s",PCFieldSplitSchurFactTypes[jac->schurfactorization]);
255: PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
256: y -= h;
257: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_USER && !jac->schur_user) {
258: PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[PC_FIELDSPLIT_SCHUR_PRE_A11]);
259: } else {
260: PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[jac->schurpre]);
261: }
262: PetscDrawStringBoxed(draw,x+wd*(cnt-1)/2.0,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
263: y -= h;
264: x = x - wd*(cnt-1)/2.0;
266: PetscDrawPushCurrentPoint(draw,x,y);
267: KSPView(jac->head->ksp,viewer);
268: PetscDrawPopCurrentPoint(draw);
269: if (jac->kspupper != jac->head->ksp) {
270: x += wd;
271: PetscDrawPushCurrentPoint(draw,x,y);
272: KSPView(jac->kspupper,viewer);
273: PetscDrawPopCurrentPoint(draw);
274: }
275: x += wd;
276: PetscDrawPushCurrentPoint(draw,x,y);
277: KSPView(jac->kspschur,viewer);
278: PetscDrawPopCurrentPoint(draw);
279: }
280: return(0);
281: }
283: static PetscErrorCode PCView_FieldSplit_GKB(PC pc,PetscViewer viewer)
284: {
285: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
286: PetscErrorCode ierr;
287: PetscBool iascii,isdraw;
288: PetscInt i,j;
289: 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: {
357: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
358: PetscInt i,nfields,*ifields,nfields_col,*ifields_col;
359: PetscBool flg,flg_col;
360: char optionname[128],splitname[8],optionname_col[128];
363: PetscMalloc1(jac->bs,&ifields);
364: PetscMalloc1(jac->bs,&ifields_col);
365: for (i=0,flg=PETSC_TRUE;; i++) {
366: PetscSNPrintf(splitname,sizeof(splitname),"%D",i);
367: PetscSNPrintf(optionname,sizeof(optionname),"-pc_fieldsplit_%D_fields",i);
368: PetscSNPrintf(optionname_col,sizeof(optionname_col),"-pc_fieldsplit_%D_fields_col",i);
369: nfields = jac->bs;
370: nfields_col = jac->bs;
371: PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,optionname,ifields,&nfields,&flg);
372: PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,optionname_col,ifields_col,&nfields_col,&flg_col);
373: if (!flg) break;
374: else if (flg && !flg_col) {
375: if (!nfields) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields");
376: PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields);
377: } else {
378: if (!nfields || !nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields");
379: if (nfields != nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Number of row and column fields must match");
380: PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields_col);
381: }
382: }
383: if (i > 0) {
384: /* Makes command-line setting of splits take precedence over setting them in code.
385: Otherwise subsequent calls to PCFieldSplitSetIS() or PCFieldSplitSetFields() would
386: create new splits, which would probably not be what the user wanted. */
387: jac->splitdefined = PETSC_TRUE;
388: }
389: PetscFree(ifields);
390: PetscFree(ifields_col);
391: return(0);
392: }
394: static PetscErrorCode PCFieldSplitSetDefaults(PC pc)
395: {
396: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
397: PetscErrorCode ierr;
398: PC_FieldSplitLink ilink = jac->head;
399: PetscBool fieldsplit_default = PETSC_FALSE,coupling = PETSC_FALSE;
400: PetscInt i;
403: /*
404: Kinda messy, but at least this now uses DMCreateFieldDecomposition().
405: Should probably be rewritten.
406: */
407: if (!ilink) {
408: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_coupling",&coupling,NULL);
409: if (pc->dm && jac->dm_splits && !jac->detect && !coupling) {
410: PetscInt numFields, f, i, j;
411: char **fieldNames;
412: IS *fields;
413: DM *dms;
414: DM subdm[128];
415: PetscBool flg;
417: DMCreateFieldDecomposition(pc->dm, &numFields, &fieldNames, &fields, &dms);
418: /* Allow the user to prescribe the splits */
419: for (i = 0, flg = PETSC_TRUE;; i++) {
420: PetscInt ifields[128];
421: IS compField;
422: char optionname[128], splitname[8];
423: PetscInt nfields = numFields;
425: PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%D_fields", i);
426: PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix, optionname, ifields, &nfields, &flg);
427: if (!flg) break;
428: if (numFields > 128) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot currently support %d > 128 fields", numFields);
429: DMCreateSubDM(pc->dm, nfields, ifields, &compField, &subdm[i]);
430: if (nfields == 1) {
431: PCFieldSplitSetIS(pc, fieldNames[ifields[0]], compField);
432: } else {
433: PetscSNPrintf(splitname, sizeof(splitname), "%D", i);
434: PCFieldSplitSetIS(pc, splitname, compField);
435: }
436: ISDestroy(&compField);
437: for (j = 0; j < nfields; ++j) {
438: f = ifields[j];
439: PetscFree(fieldNames[f]);
440: ISDestroy(&fields[f]);
441: }
442: }
443: if (i == 0) {
444: for (f = 0; f < numFields; ++f) {
445: PCFieldSplitSetIS(pc, fieldNames[f], fields[f]);
446: PetscFree(fieldNames[f]);
447: ISDestroy(&fields[f]);
448: }
449: } else {
450: for (j=0; j<numFields; j++) {
451: DMDestroy(dms+j);
452: }
453: PetscFree(dms);
454: PetscMalloc1(i, &dms);
455: for (j = 0; j < i; ++j) dms[j] = subdm[j];
456: }
457: PetscFree(fieldNames);
458: PetscFree(fields);
459: if (dms) {
460: PetscInfo(pc, "Setting up physics based fieldsplit preconditioner using the embedded DM\n");
461: for (ilink = jac->head, i = 0; ilink; ilink = ilink->next, ++i) {
462: const char *prefix;
463: PetscObjectGetOptionsPrefix((PetscObject)(ilink->ksp),&prefix);
464: PetscObjectSetOptionsPrefix((PetscObject)(dms[i]), prefix);
465: KSPSetDM(ilink->ksp, dms[i]);
466: KSPSetDMActive(ilink->ksp, PETSC_FALSE);
467: {
468: PetscErrorCode (*func)(KSP,Mat,Mat,void*);
469: void *ctx;
471: DMKSPGetComputeOperators(pc->dm, &func, &ctx);
472: DMKSPSetComputeOperators(dms[i], func, ctx);
473: }
474: PetscObjectIncrementTabLevel((PetscObject)dms[i],(PetscObject)ilink->ksp,0);
475: DMDestroy(&dms[i]);
476: }
477: PetscFree(dms);
478: }
479: } else {
480: if (jac->bs <= 0) {
481: if (pc->pmat) {
482: MatGetBlockSize(pc->pmat,&jac->bs);
483: } else jac->bs = 1;
484: }
486: if (jac->detect) {
487: IS zerodiags,rest;
488: PetscInt nmin,nmax;
490: MatGetOwnershipRange(pc->mat,&nmin,&nmax);
491: if (jac->diag_use_amat) {
492: MatFindZeroDiagonals(pc->mat,&zerodiags);
493: } else {
494: MatFindZeroDiagonals(pc->pmat,&zerodiags);
495: }
496: ISComplement(zerodiags,nmin,nmax,&rest);
497: PCFieldSplitSetIS(pc,"0",rest);
498: PCFieldSplitSetIS(pc,"1",zerodiags);
499: ISDestroy(&zerodiags);
500: ISDestroy(&rest);
501: } else if (coupling) {
502: IS coupling,rest;
503: PetscInt nmin,nmax;
505: MatGetOwnershipRange(pc->mat,&nmin,&nmax);
506: if (jac->offdiag_use_amat) {
507: MatFindOffBlockDiagonalEntries(pc->mat,&coupling);
508: } else {
509: MatFindOffBlockDiagonalEntries(pc->pmat,&coupling);
510: }
511: ISCreateStride(PetscObjectComm((PetscObject)pc->mat),nmax-nmin,nmin,1,&rest);
512: ISSetIdentity(rest);
513: PCFieldSplitSetIS(pc,"0",rest);
514: PCFieldSplitSetIS(pc,"1",coupling);
515: ISDestroy(&coupling);
516: ISDestroy(&rest);
517: } else {
518: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&fieldsplit_default,NULL);
519: if (!fieldsplit_default) {
520: /* Allow user to set fields from command line, if bs was known at the time of PCSetFromOptions_FieldSplit()
521: then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */
522: PCFieldSplitSetRuntimeSplits_Private(pc);
523: if (jac->splitdefined) {PetscInfo(pc,"Splits defined using the options database\n");}
524: }
525: if ((fieldsplit_default || !jac->splitdefined) && !jac->isrestrict) {
526: Mat M = pc->pmat;
527: PetscBool isnest;
529: PetscInfo(pc,"Using default splitting of fields\n");
530: PetscObjectTypeCompare((PetscObject)pc->pmat,MATNEST,&isnest);
531: if (!isnest) {
532: M = pc->mat;
533: PetscObjectTypeCompare((PetscObject)pc->mat,MATNEST,&isnest);
534: }
535: if (isnest) {
536: IS *fields;
537: PetscInt nf;
539: MatNestGetSize(M,&nf,NULL);
540: PetscMalloc1(nf,&fields);
541: MatNestGetISs(M,fields,NULL);
542: for (i=0;i<nf;i++) {
543: PCFieldSplitSetIS(pc,NULL,fields[i]);
544: }
545: PetscFree(fields);
546: } else {
547: for (i=0; i<jac->bs; i++) {
548: char splitname[8];
549: PetscSNPrintf(splitname,sizeof(splitname),"%D",i);
550: PCFieldSplitSetFields(pc,splitname,1,&i,&i);
551: }
552: jac->defaultsplit = PETSC_TRUE;
553: }
554: }
555: }
556: }
557: } else if (jac->nsplits == 1) {
558: if (ilink->is) {
559: IS is2;
560: PetscInt nmin,nmax;
562: MatGetOwnershipRange(pc->mat,&nmin,&nmax);
563: ISComplement(ilink->is,nmin,nmax,&is2);
564: PCFieldSplitSetIS(pc,"1",is2);
565: ISDestroy(&is2);
566: } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Must provide at least two sets of fields to PCFieldSplit()");
567: }
569: if (jac->nsplits < 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Unhandled case, must have at least two fields, not %d", jac->nsplits);
570: return(0);
571: }
573: static PetscErrorCode MatGolubKahanComputeExplicitOperator(Mat A,Mat B,Mat C,Mat *H,PetscReal gkbnu)
574: {
575: PetscErrorCode ierr;
576: Mat BT,T;
577: PetscReal nrmT,nrmB;
580: MatHermitianTranspose(C,MAT_INITIAL_MATRIX,&T); /* Test if augmented matrix is symmetric */
581: MatAXPY(T,-1.0,B,DIFFERENT_NONZERO_PATTERN);
582: MatNorm(T,NORM_1,&nrmT);
583: MatNorm(B,NORM_1,&nrmB);
584: if (nrmB > 0) {
585: if (nrmT/nrmB >= PETSC_SMALL) {
586: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix is not symmetric/hermitian, GKB is not applicable.");
587: }
588: }
589: /* Compute augmented Lagrangian matrix H = A00 + nu*A01*A01'. This corresponds to */
590: /* setting N := 1/nu*I in [Ar13]. */
591: MatHermitianTranspose(B,MAT_INITIAL_MATRIX,&BT);
592: MatMatMult(B,BT,MAT_INITIAL_MATRIX,PETSC_DEFAULT,H); /* H = A01*A01' */
593: MatAYPX(*H,gkbnu,A,DIFFERENT_NONZERO_PATTERN); /* H = A00 + nu*A01*A01' */
595: MatDestroy(&BT);
596: MatDestroy(&T);
597: return(0);
598: }
600: PETSC_EXTERN PetscErrorCode PetscOptionsFindPairPrefix_Private(PetscOptions,const char pre[], const char name[],const char *value[],PetscBool *flg);
602: static PetscErrorCode PCSetUp_FieldSplit(PC pc)
603: {
604: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
605: PetscErrorCode ierr;
606: PC_FieldSplitLink ilink;
607: PetscInt i,nsplit;
608: PetscBool sorted, sorted_col;
611: pc->failedreason = PC_NOERROR;
612: PCFieldSplitSetDefaults(pc);
613: nsplit = jac->nsplits;
614: ilink = jac->head;
616: /* get the matrices for each split */
617: if (!jac->issetup) {
618: PetscInt rstart,rend,nslots,bs;
620: jac->issetup = PETSC_TRUE;
622: /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */
623: if (jac->defaultsplit || !ilink->is) {
624: if (jac->bs <= 0) jac->bs = nsplit;
625: }
626: bs = jac->bs;
627: MatGetOwnershipRange(pc->pmat,&rstart,&rend);
628: nslots = (rend - rstart)/bs;
629: for (i=0; i<nsplit; i++) {
630: if (jac->defaultsplit) {
631: ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+i,nsplit,&ilink->is);
632: ISDuplicate(ilink->is,&ilink->is_col);
633: } else if (!ilink->is) {
634: if (ilink->nfields > 1) {
635: PetscInt *ii,*jj,j,k,nfields = ilink->nfields,*fields = ilink->fields,*fields_col = ilink->fields_col;
636: PetscMalloc1(ilink->nfields*nslots,&ii);
637: PetscMalloc1(ilink->nfields*nslots,&jj);
638: for (j=0; j<nslots; j++) {
639: for (k=0; k<nfields; k++) {
640: ii[nfields*j + k] = rstart + bs*j + fields[k];
641: jj[nfields*j + k] = rstart + bs*j + fields_col[k];
642: }
643: }
644: ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,ii,PETSC_OWN_POINTER,&ilink->is);
645: ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,jj,PETSC_OWN_POINTER,&ilink->is_col);
646: ISSetBlockSize(ilink->is, nfields);
647: ISSetBlockSize(ilink->is_col, nfields);
648: } else {
649: ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields[0],bs,&ilink->is);
650: ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields_col[0],bs,&ilink->is_col);
651: }
652: }
653: ISSorted(ilink->is,&sorted);
654: if (ilink->is_col) { ISSorted(ilink->is_col,&sorted_col); }
655: if (!sorted || !sorted_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Fields must be sorted when creating split");
656: ilink = ilink->next;
657: }
658: }
660: ilink = jac->head;
661: if (!jac->pmat) {
662: Vec xtmp;
664: MatCreateVecs(pc->pmat,&xtmp,NULL);
665: PetscMalloc1(nsplit,&jac->pmat);
666: PetscMalloc2(nsplit,&jac->x,nsplit,&jac->y);
667: for (i=0; i<nsplit; i++) {
668: MatNullSpace sp;
670: /* Check for preconditioning matrix attached to IS */
671: PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &jac->pmat[i]);
672: if (jac->pmat[i]) {
673: PetscObjectReference((PetscObject) jac->pmat[i]);
674: if (jac->type == PC_COMPOSITE_SCHUR) {
675: jac->schur_user = jac->pmat[i];
677: PetscObjectReference((PetscObject) jac->schur_user);
678: }
679: } else {
680: const char *prefix;
681: MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->pmat[i]);
682: KSPGetOptionsPrefix(ilink->ksp,&prefix);
683: MatSetOptionsPrefix(jac->pmat[i],prefix);
684: MatViewFromOptions(jac->pmat[i],NULL,"-mat_view");
685: }
686: /* create work vectors for each split */
687: MatCreateVecs(jac->pmat[i],&jac->x[i],&jac->y[i]);
688: ilink->x = jac->x[i]; ilink->y = jac->y[i]; ilink->z = NULL;
689: /* compute scatter contexts needed by multiplicative versions and non-default splits */
690: VecScatterCreate(xtmp,ilink->is,jac->x[i],NULL,&ilink->sctx);
691: PetscObjectQuery((PetscObject) ilink->is, "nearnullspace", (PetscObject*) &sp);
692: if (sp) {
693: MatSetNearNullSpace(jac->pmat[i], sp);
694: }
695: ilink = ilink->next;
696: }
697: VecDestroy(&xtmp);
698: } else {
699: MatReuse scall;
700: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
701: for (i=0; i<nsplit; i++) {
702: MatDestroy(&jac->pmat[i]);
703: }
704: scall = MAT_INITIAL_MATRIX;
705: } else scall = MAT_REUSE_MATRIX;
707: for (i=0; i<nsplit; i++) {
708: Mat pmat;
710: /* Check for preconditioning matrix attached to IS */
711: PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &pmat);
712: if (!pmat) {
713: MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is_col,scall,&jac->pmat[i]);
714: }
715: ilink = ilink->next;
716: }
717: }
718: if (jac->diag_use_amat) {
719: ilink = jac->head;
720: if (!jac->mat) {
721: PetscMalloc1(nsplit,&jac->mat);
722: for (i=0; i<nsplit; i++) {
723: MatCreateSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->mat[i]);
724: ilink = ilink->next;
725: }
726: } else {
727: MatReuse scall;
728: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
729: for (i=0; i<nsplit; i++) {
730: MatDestroy(&jac->mat[i]);
731: }
732: scall = MAT_INITIAL_MATRIX;
733: } else scall = MAT_REUSE_MATRIX;
735: for (i=0; i<nsplit; i++) {
736: MatCreateSubMatrix(pc->mat,ilink->is,ilink->is_col,scall,&jac->mat[i]);
737: ilink = ilink->next;
738: }
739: }
740: } else {
741: jac->mat = jac->pmat;
742: }
744: /* Check for null space attached to IS */
745: ilink = jac->head;
746: for (i=0; i<nsplit; i++) {
747: MatNullSpace sp;
749: PetscObjectQuery((PetscObject) ilink->is, "nullspace", (PetscObject*) &sp);
750: if (sp) {
751: MatSetNullSpace(jac->mat[i], sp);
752: }
753: ilink = ilink->next;
754: }
756: if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR && jac->type != PC_COMPOSITE_GKB) {
757: /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */
758: /* FIXME: Can/should we reuse jac->mat whenever (jac->diag_use_amat) is true? */
759: ilink = jac->head;
760: if (nsplit == 2 && jac->type == PC_COMPOSITE_MULTIPLICATIVE) {
761: /* 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 */
762: if (!jac->Afield) {
763: PetscCalloc1(nsplit,&jac->Afield);
764: if (jac->offdiag_use_amat) {
765: MatCreateSubMatrix(pc->mat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);
766: } else {
767: MatCreateSubMatrix(pc->pmat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);
768: }
769: } else {
770: MatReuse scall;
772: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
773: MatDestroy(&jac->Afield[1]);
774: scall = MAT_INITIAL_MATRIX;
775: } else scall = MAT_REUSE_MATRIX;
777: if (jac->offdiag_use_amat) {
778: MatCreateSubMatrix(pc->mat,ilink->next->is,ilink->is,scall,&jac->Afield[1]);
779: } else {
780: MatCreateSubMatrix(pc->pmat,ilink->next->is,ilink->is,scall,&jac->Afield[1]);
781: }
782: }
783: } else {
784: if (!jac->Afield) {
785: PetscMalloc1(nsplit,&jac->Afield);
786: for (i=0; i<nsplit; i++) {
787: if (jac->offdiag_use_amat) {
788: MatCreateSubMatrix(pc->mat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);
789: } else {
790: MatCreateSubMatrix(pc->pmat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);
791: }
792: ilink = ilink->next;
793: }
794: } else {
795: MatReuse scall;
796: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
797: for (i=0; i<nsplit; i++) {
798: MatDestroy(&jac->Afield[i]);
799: }
800: scall = MAT_INITIAL_MATRIX;
801: } else scall = MAT_REUSE_MATRIX;
803: for (i=0; i<nsplit; i++) {
804: if (jac->offdiag_use_amat) {
805: MatCreateSubMatrix(pc->mat,ilink->is,NULL,scall,&jac->Afield[i]);
806: } else {
807: MatCreateSubMatrix(pc->pmat,ilink->is,NULL,scall,&jac->Afield[i]);
808: }
809: ilink = ilink->next;
810: }
811: }
812: }
813: }
815: if (jac->type == PC_COMPOSITE_SCHUR) {
816: IS ccis;
817: PetscBool isspd;
818: PetscInt rstart,rend;
819: char lscname[256];
820: PetscObject LSC_L;
822: if (nsplit != 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"To use Schur complement preconditioner you must have exactly 2 fields");
824: /* If pc->mat is SPD, don't scale by -1 the Schur complement */
825: if (jac->schurscale == (PetscScalar)-1.0) {
826: MatGetOption(pc->pmat,MAT_SPD,&isspd);
827: jac->schurscale = (isspd == PETSC_TRUE) ? 1.0 : -1.0;
828: }
830: /* When extracting off-diagonal submatrices, we take complements from this range */
831: MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);
833: if (jac->schur) {
834: KSP kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper;
835: MatReuse scall;
837: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
838: scall = MAT_INITIAL_MATRIX;
839: MatDestroy(&jac->B);
840: MatDestroy(&jac->C);
841: } else scall = MAT_REUSE_MATRIX;
843: MatSchurComplementGetKSP(jac->schur, &kspInner);
844: ilink = jac->head;
845: ISComplement(ilink->is_col,rstart,rend,&ccis);
846: if (jac->offdiag_use_amat) {
847: MatCreateSubMatrix(pc->mat,ilink->is,ccis,scall,&jac->B);
848: } else {
849: MatCreateSubMatrix(pc->pmat,ilink->is,ccis,scall,&jac->B);
850: }
851: ISDestroy(&ccis);
852: ilink = ilink->next;
853: ISComplement(ilink->is_col,rstart,rend,&ccis);
854: if (jac->offdiag_use_amat) {
855: MatCreateSubMatrix(pc->mat,ilink->is,ccis,scall,&jac->C);
856: } else {
857: MatCreateSubMatrix(pc->pmat,ilink->is,ccis,scall,&jac->C);
858: }
859: ISDestroy(&ccis);
860: MatSchurComplementUpdateSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);
861: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
862: MatDestroy(&jac->schurp);
863: MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);
864: }
865: if (kspA != kspInner) {
866: KSPSetOperators(kspA,jac->mat[0],jac->pmat[0]);
867: }
868: if (kspUpper != kspA) {
869: KSPSetOperators(kspUpper,jac->mat[0],jac->pmat[0]);
870: }
871: KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));
872: } else {
873: const char *Dprefix;
874: char schurprefix[256], schurmatprefix[256];
875: char schurtestoption[256];
876: MatNullSpace sp;
877: PetscBool flg;
878: KSP kspt;
880: /* extract the A01 and A10 matrices */
881: ilink = jac->head;
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->B);
885: } else {
886: MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
887: }
888: ISDestroy(&ccis);
889: ilink = ilink->next;
890: ISComplement(ilink->is_col,rstart,rend,&ccis);
891: if (jac->offdiag_use_amat) {
892: MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
893: } else {
894: MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
895: }
896: ISDestroy(&ccis);
898: /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */
899: MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur);
900: MatSetType(jac->schur,MATSCHURCOMPLEMENT);
901: MatSchurComplementSetSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);
902: PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
903: MatSetOptionsPrefix(jac->schur,schurmatprefix);
904: MatSchurComplementGetKSP(jac->schur,&kspt);
905: KSPSetOptionsPrefix(kspt,schurmatprefix);
907: /* Note: this is not true in general */
908: MatGetNullSpace(jac->mat[1], &sp);
909: if (sp) {
910: MatSetNullSpace(jac->schur, sp);
911: }
913: PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);
914: PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
915: if (flg) {
916: DM dmInner;
917: KSP kspInner;
918: PC pcInner;
920: MatSchurComplementGetKSP(jac->schur, &kspInner);
921: KSPReset(kspInner);
922: KSPSetOperators(kspInner,jac->mat[0],jac->pmat[0]);
923: PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
924: /* Indent this deeper to emphasize the "inner" nature of this solver. */
925: PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject) pc, 2);
926: PetscObjectIncrementTabLevel((PetscObject)kspInner->pc, (PetscObject) pc, 2);
927: KSPSetOptionsPrefix(kspInner, schurprefix);
929: /* Set DM for new solver */
930: KSPGetDM(jac->head->ksp, &dmInner);
931: KSPSetDM(kspInner, dmInner);
932: KSPSetDMActive(kspInner, PETSC_FALSE);
934: /* Defaults to PCKSP as preconditioner */
935: KSPGetPC(kspInner, &pcInner);
936: PCSetType(pcInner, PCKSP);
937: PCKSPSetKSP(pcInner, jac->head->ksp);
938: } else {
939: /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or
940: * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact,
941: * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for
942: * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make
943: * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used
944: * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */
945: KSPSetType(jac->head->ksp,KSPGMRES);
946: MatSchurComplementSetKSP(jac->schur,jac->head->ksp);
947: }
948: KSPSetOperators(jac->head->ksp,jac->mat[0],jac->pmat[0]);
949: KSPSetFromOptions(jac->head->ksp);
950: MatSetFromOptions(jac->schur);
952: PetscObjectTypeCompare((PetscObject)jac->schur, MATSCHURCOMPLEMENT, &flg);
953: if (flg) { /* Need to do this otherwise PCSetUp_KSP will overwrite the amat of jac->head->ksp */
954: KSP kspInner;
955: PC pcInner;
957: MatSchurComplementGetKSP(jac->schur, &kspInner);
958: KSPGetPC(kspInner, &pcInner);
959: PetscObjectTypeCompare((PetscObject)pcInner, PCKSP, &flg);
960: if (flg) {
961: KSP ksp;
963: PCKSPGetKSP(pcInner, &ksp);
964: if (ksp == jac->head->ksp) {
965: PCSetUseAmat(pcInner, PETSC_TRUE);
966: }
967: }
968: }
969: PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);
970: PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
971: if (flg) {
972: DM dmInner;
974: PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
975: KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper);
976: KSPSetErrorIfNotConverged(jac->kspupper,pc->erroriffailure);
977: KSPSetOptionsPrefix(jac->kspupper, schurprefix);
978: PetscObjectIncrementTabLevel((PetscObject)jac->kspupper, (PetscObject) pc, 1);
979: PetscObjectIncrementTabLevel((PetscObject)jac->kspupper->pc, (PetscObject) pc, 1);
980: KSPGetDM(jac->head->ksp, &dmInner);
981: KSPSetDM(jac->kspupper, dmInner);
982: KSPSetDMActive(jac->kspupper, PETSC_FALSE);
983: KSPSetFromOptions(jac->kspupper);
984: KSPSetOperators(jac->kspupper,jac->mat[0],jac->pmat[0]);
985: VecDuplicate(jac->head->x, &jac->head->z);
986: } else {
987: jac->kspupper = jac->head->ksp;
988: PetscObjectReference((PetscObject) jac->head->ksp);
989: }
991: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
992: MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);
993: }
994: KSPCreate(PetscObjectComm((PetscObject)pc),&jac->kspschur);
995: KSPSetErrorIfNotConverged(jac->kspschur,pc->erroriffailure);
996: PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);
997: PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);
998: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) {
999: PC pcschur;
1000: KSPGetPC(jac->kspschur,&pcschur);
1001: PCSetType(pcschur,PCNONE);
1002: /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */
1003: } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) {
1004: MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user);
1005: }
1006: KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));
1007: KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix);
1008: KSPSetOptionsPrefix(jac->kspschur, Dprefix);
1009: /* propagate DM */
1010: {
1011: DM sdm;
1012: KSPGetDM(jac->head->next->ksp, &sdm);
1013: if (sdm) {
1014: KSPSetDM(jac->kspschur, sdm);
1015: KSPSetDMActive(jac->kspschur, PETSC_FALSE);
1016: }
1017: }
1018: /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
1019: /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */
1020: KSPSetFromOptions(jac->kspschur);
1021: }
1022: MatAssemblyBegin(jac->schur,MAT_FINAL_ASSEMBLY);
1023: MatAssemblyEnd(jac->schur,MAT_FINAL_ASSEMBLY);
1025: /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */
1026: PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_L",ilink->splitname);
1027: PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);
1028: if (!LSC_L) {PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);}
1029: if (LSC_L) {PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L);}
1030: PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_Lp",ilink->splitname);
1031: PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);
1032: if (!LSC_L) {PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);}
1033: if (LSC_L) {PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L);}
1034: } else if (jac->type == PC_COMPOSITE_GKB) {
1035: IS ccis;
1036: PetscInt rstart,rend;
1038: if (nsplit != 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"To use GKB preconditioner you must have exactly 2 fields");
1040: ilink = jac->head;
1042: /* When extracting off-diagonal submatrices, we take complements from this range */
1043: MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);
1045: ISComplement(ilink->is_col,rstart,rend,&ccis);
1046: if (jac->offdiag_use_amat) {
1047: MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
1048: } else {
1049: MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
1050: }
1051: ISDestroy(&ccis);
1052: /* Create work vectors for GKB algorithm */
1053: VecDuplicate(ilink->x,&jac->u);
1054: VecDuplicate(ilink->x,&jac->Hu);
1055: VecDuplicate(ilink->x,&jac->w2);
1056: ilink = ilink->next;
1057: ISComplement(ilink->is_col,rstart,rend,&ccis);
1058: if (jac->offdiag_use_amat) {
1059: MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
1060: } else {
1061: MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
1062: }
1063: ISDestroy(&ccis);
1064: /* Create work vectors for GKB algorithm */
1065: VecDuplicate(ilink->x,&jac->v);
1066: VecDuplicate(ilink->x,&jac->d);
1067: VecDuplicate(ilink->x,&jac->w1);
1068: MatGolubKahanComputeExplicitOperator(jac->mat[0],jac->B,jac->C,&jac->H,jac->gkbnu);
1069: PetscCalloc1(jac->gkbdelay,&jac->vecz);
1071: ilink = jac->head;
1072: KSPSetOperators(ilink->ksp,jac->H,jac->H);
1073: if (!jac->suboptionsset) {KSPSetFromOptions(ilink->ksp);}
1074: /* Create gkb_monitor context */
1075: if (jac->gkbmonitor) {
1076: PetscInt tablevel;
1077: PetscViewerCreate(PETSC_COMM_WORLD,&jac->gkbviewer);
1078: PetscViewerSetType(jac->gkbviewer,PETSCVIEWERASCII);
1079: PetscObjectGetTabLevel((PetscObject)ilink->ksp,&tablevel);
1080: PetscViewerASCIISetTab(jac->gkbviewer,tablevel);
1081: PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)ilink->ksp,1);
1082: }
1083: } else {
1084: /* set up the individual splits' PCs */
1085: i = 0;
1086: ilink = jac->head;
1087: while (ilink) {
1088: KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i]);
1089: /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
1090: if (!jac->suboptionsset) {KSPSetFromOptions(ilink->ksp);}
1091: i++;
1092: ilink = ilink->next;
1093: }
1094: }
1096: jac->suboptionsset = PETSC_TRUE;
1097: return(0);
1098: }
1100: #define FieldSplitSplitSolveAdd(ilink,xx,yy) \
1101: (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
1102: VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
1103: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\
1104: KSPSolve(ilink->ksp,ilink->x,ilink->y) || \
1105: KSPCheckSolve(ilink->ksp,pc,ilink->y) || \
1106: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\
1107: VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \
1108: VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE))
1110: static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y)
1111: {
1112: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1113: PetscErrorCode ierr;
1114: PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
1115: KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;
1118: switch (jac->schurfactorization) {
1119: case PC_FIELDSPLIT_SCHUR_FACT_DIAG:
1120: /* [A00 0; 0 -S], positive definite, suitable for MINRES */
1121: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1122: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
1123: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1124: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1125: KSPSolve(kspA,ilinkA->x,ilinkA->y);
1126: KSPCheckSolve(kspA,pc,ilinkA->y);
1127: PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1128: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1129: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
1130: PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1131: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
1132: KSPCheckSolve(jac->kspschur,pc,ilinkD->y);
1133: PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1134: VecScale(ilinkD->y,jac->schurscale);
1135: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1136: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1137: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1138: break;
1139: case PC_FIELDSPLIT_SCHUR_FACT_LOWER:
1140: /* [A00 0; A10 S], suitable for left preconditioning */
1141: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1142: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1143: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1144: KSPSolve(kspA,ilinkA->x,ilinkA->y);
1145: KSPCheckSolve(kspA,pc,ilinkA->y);
1146: PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1147: MatMult(jac->C,ilinkA->y,ilinkD->x);
1148: VecScale(ilinkD->x,-1.);
1149: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
1150: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1151: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
1152: PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1153: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
1154: KSPCheckSolve(jac->kspschur,pc,ilinkD->y);
1155: PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1156: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1157: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1158: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1159: break;
1160: case PC_FIELDSPLIT_SCHUR_FACT_UPPER:
1161: /* [A00 A01; 0 S], suitable for right preconditioning */
1162: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
1163: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_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); MatMult(jac->B,ilinkD->y,ilinkA->x);
1168: VecScale(ilinkA->x,-1.);
1169: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);
1170: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1171: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);
1172: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1173: KSPSolve(kspA,ilinkA->x,ilinkA->y);
1174: KSPCheckSolve(kspA,pc,ilinkA->y);
1175: PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1176: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1177: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1178: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1179: break;
1180: case PC_FIELDSPLIT_SCHUR_FACT_FULL:
1181: /* [1 0; A10 A00^{-1} 1] [A00 0; 0 S] [1 A00^{-1}A01; 0 1] */
1182: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1183: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1184: PetscLogEventBegin(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);
1185: KSPSolve(kspLower,ilinkA->x,ilinkA->y);
1186: KSPCheckSolve(kspLower,pc,ilinkA->y);
1187: PetscLogEventEnd(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);
1188: MatMult(jac->C,ilinkA->y,ilinkD->x);
1189: VecScale(ilinkD->x,-1.0);
1190: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
1191: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
1193: PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1194: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
1195: KSPCheckSolve(jac->kspschur,pc,ilinkD->y);
1196: PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1197: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1199: if (kspUpper == kspA) {
1200: MatMult(jac->B,ilinkD->y,ilinkA->y);
1201: VecAXPY(ilinkA->x,-1.0,ilinkA->y);
1202: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1203: KSPSolve(kspA,ilinkA->x,ilinkA->y);
1204: KSPCheckSolve(kspA,pc,ilinkA->y);
1205: PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1206: } else {
1207: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1208: KSPSolve(kspA,ilinkA->x,ilinkA->y);
1209: KSPCheckSolve(kspA,pc,ilinkA->y);
1210: MatMult(jac->B,ilinkD->y,ilinkA->x);
1211: PetscLogEventBegin(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);
1212: KSPSolve(kspUpper,ilinkA->x,ilinkA->z);
1213: KSPCheckSolve(kspUpper,pc,ilinkA->z);
1214: PetscLogEventEnd(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);
1215: VecAXPY(ilinkA->y,-1.0,ilinkA->z);
1216: }
1217: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1218: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1219: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1220: }
1221: return(0);
1222: }
1224: static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y)
1225: {
1226: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1227: PetscErrorCode ierr;
1228: PC_FieldSplitLink ilink = jac->head;
1229: PetscInt cnt,bs;
1232: if (jac->type == PC_COMPOSITE_ADDITIVE) {
1233: if (jac->defaultsplit) {
1234: VecGetBlockSize(x,&bs);
1235: if (jac->bs > 0 && bs != jac->bs) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Blocksize of x vector %D does not match fieldsplit blocksize %D",bs,jac->bs);
1236: VecGetBlockSize(y,&bs);
1237: if (jac->bs > 0 && bs != jac->bs) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Blocksize of y vector %D does not match fieldsplit blocksize %D",bs,jac->bs);
1238: VecStrideGatherAll(x,jac->x,INSERT_VALUES);
1239: while (ilink) {
1240: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1241: KSPSolve(ilink->ksp,ilink->x,ilink->y);
1242: KSPCheckSolve(ilink->ksp,pc,ilink->y);
1243: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1244: ilink = ilink->next;
1245: }
1246: VecStrideScatterAll(jac->y,y,INSERT_VALUES);
1247: } else {
1248: VecSet(y,0.0);
1249: while (ilink) {
1250: FieldSplitSplitSolveAdd(ilink,x,y);
1251: ilink = ilink->next;
1252: }
1253: }
1254: } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) {
1255: VecSet(y,0.0);
1256: /* solve on first block for first block variables */
1257: VecScatterBegin(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);
1258: VecScatterEnd(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);
1259: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1260: KSPSolve(ilink->ksp,ilink->x,ilink->y);
1261: KSPCheckSolve(ilink->ksp,pc,ilink->y);
1262: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1263: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1264: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1266: /* compute the residual only onto second block variables using first block variables */
1267: MatMult(jac->Afield[1],ilink->y,ilink->next->x);
1268: ilink = ilink->next;
1269: VecScale(ilink->x,-1.0);
1270: VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1271: VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1273: /* solve on second block variables */
1274: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1275: KSPSolve(ilink->ksp,ilink->x,ilink->y);
1276: KSPCheckSolve(ilink->ksp,pc,ilink->y);
1277: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1278: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1279: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1280: } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1281: if (!jac->w1) {
1282: VecDuplicate(x,&jac->w1);
1283: VecDuplicate(x,&jac->w2);
1284: }
1285: VecSet(y,0.0);
1286: FieldSplitSplitSolveAdd(ilink,x,y);
1287: cnt = 1;
1288: while (ilink->next) {
1289: ilink = ilink->next;
1290: /* compute the residual only over the part of the vector needed */
1291: MatMult(jac->Afield[cnt++],y,ilink->x);
1292: VecScale(ilink->x,-1.0);
1293: VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1294: VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1295: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1296: KSPSolve(ilink->ksp,ilink->x,ilink->y);
1297: KSPCheckSolve(ilink->ksp,pc,ilink->y);
1298: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1299: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1300: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1301: }
1302: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1303: cnt -= 2;
1304: while (ilink->previous) {
1305: ilink = ilink->previous;
1306: /* compute the residual only over the part of the vector needed */
1307: MatMult(jac->Afield[cnt--],y,ilink->x);
1308: VecScale(ilink->x,-1.0);
1309: VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1310: VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1311: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1312: KSPSolve(ilink->ksp,ilink->x,ilink->y);
1313: KSPCheckSolve(ilink->ksp,pc,ilink->y);
1314: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1315: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1316: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1317: }
1318: }
1319: } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type);
1320: return(0);
1321: }
1323: static PetscErrorCode PCApply_FieldSplit_GKB(PC pc,Vec x,Vec y)
1324: {
1325: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1326: PetscErrorCode ierr;
1327: PC_FieldSplitLink ilinkA = jac->head,ilinkD = ilinkA->next;
1328: KSP ksp = ilinkA->ksp;
1329: Vec u,v,Hu,d,work1,work2;
1330: PetscScalar alpha,z,nrmz2,*vecz;
1331: PetscReal lowbnd,nu,beta;
1332: PetscInt j,iterGKB;
1335: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1336: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
1337: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1338: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
1340: u = jac->u;
1341: v = jac->v;
1342: Hu = jac->Hu;
1343: d = jac->d;
1344: work1 = jac->w1;
1345: work2 = jac->w2;
1346: vecz = jac->vecz;
1348: /* Change RHS to comply with matrix regularization H = A + nu*B*B' */
1349: /* Add q = q + nu*B*b */
1350: if (jac->gkbnu) {
1351: nu = jac->gkbnu;
1352: VecScale(ilinkD->x,jac->gkbnu);
1353: MatMultAdd(jac->B,ilinkD->x,ilinkA->x,ilinkA->x); /* q = q + nu*B*b */
1354: } else {
1355: /* Situation when no augmented Lagrangian is used. Then we set inner */
1356: /* matrix N = I in [Ar13], and thus nu = 1. */
1357: nu = 1;
1358: }
1360: /* Transform rhs from [q,tilde{b}] to [0,b] */
1361: PetscLogEventBegin(ilinkA->event,ksp,ilinkA->x,ilinkA->y,NULL);
1362: KSPSolve(ksp,ilinkA->x,ilinkA->y);
1363: KSPCheckSolve(ksp,pc,ilinkA->y);
1364: PetscLogEventEnd(ilinkA->event,ksp,ilinkA->x,ilinkA->y,NULL);
1365: MatMultHermitianTranspose(jac->B,ilinkA->y,work1);
1366: VecAXPBY(work1,1.0/nu,-1.0,ilinkD->x); /* c = b - B'*x */
1368: /* First step of algorithm */
1369: VecNorm(work1,NORM_2,&beta); /* beta = sqrt(nu*c'*c)*/
1370: KSPCheckDot(ksp,beta);
1371: beta = PetscSqrtReal(nu)*beta;
1372: VecAXPBY(v,nu/beta,0.0,work1); /* v = nu/beta *c */
1373: MatMult(jac->B,v,work2); /* u = H^{-1}*B*v */
1374: PetscLogEventBegin(ilinkA->event,ksp,work2,u,NULL);
1375: KSPSolve(ksp,work2,u);
1376: KSPCheckSolve(ksp,pc,u);
1377: PetscLogEventEnd(ilinkA->event,ksp,work2,u,NULL);
1378: MatMult(jac->H,u,Hu); /* alpha = u'*H*u */
1379: VecDot(Hu,u,&alpha);
1380: KSPCheckDot(ksp,alpha);
1381: if (PetscRealPart(alpha) <= 0.0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_NOT_CONVERGED,"GKB preconditioner diverged, H is not positive definite");
1382: alpha = PetscSqrtReal(PetscAbsScalar(alpha));
1383: VecScale(u,1.0/alpha);
1384: VecAXPBY(d,1.0/alpha,0.0,v); /* v = nu/beta *c */
1386: z = beta/alpha;
1387: vecz[1] = z;
1389: /* Computation of first iterate x(1) and p(1) */
1390: VecAXPY(ilinkA->y,z,u);
1391: VecCopy(d,ilinkD->y);
1392: VecScale(ilinkD->y,-z);
1394: iterGKB = 1; lowbnd = 2*jac->gkbtol;
1395: if (jac->gkbmonitor) {
1396: PetscViewerASCIIPrintf(jac->gkbviewer,"%3D GKB Lower bound estimate %14.12e\n",iterGKB,lowbnd);
1397: }
1399: while (iterGKB < jac->gkbmaxit && lowbnd > jac->gkbtol) {
1400: iterGKB += 1;
1401: MatMultHermitianTranspose(jac->B,u,work1); /* v <- nu*(B'*u-alpha/nu*v) */
1402: VecAXPBY(v,nu,-alpha,work1);
1403: VecNorm(v,NORM_2,&beta); /* beta = sqrt(nu)*v'*v */
1404: beta = beta/PetscSqrtReal(nu);
1405: VecScale(v,1.0/beta);
1406: MatMult(jac->B,v,work2); /* u <- H^{-1}*(B*v-beta*H*u) */
1407: MatMult(jac->H,u,Hu);
1408: VecAXPY(work2,-beta,Hu);
1409: PetscLogEventBegin(ilinkA->event,ksp,work2,u,NULL);
1410: KSPSolve(ksp,work2,u);
1411: KSPCheckSolve(ksp,pc,u);
1412: PetscLogEventEnd(ilinkA->event,ksp,work2,u,NULL);
1413: MatMult(jac->H,u,Hu); /* alpha = u'*H*u */
1414: VecDot(Hu,u,&alpha);
1415: KSPCheckDot(ksp,alpha);
1416: if (PetscRealPart(alpha) <= 0.0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_NOT_CONVERGED,"GKB preconditioner diverged, H is not positive definite");
1417: alpha = PetscSqrtReal(PetscAbsScalar(alpha));
1418: VecScale(u,1.0/alpha);
1420: z = -beta/alpha*z; /* z <- beta/alpha*z */
1421: vecz[0] = z;
1423: /* Computation of new iterate x(i+1) and p(i+1) */
1424: VecAXPBY(d,1.0/alpha,-beta/alpha,v); /* d = (v-beta*d)/alpha */
1425: VecAXPY(ilinkA->y,z,u); /* r = r + z*u */
1426: VecAXPY(ilinkD->y,-z,d); /* p = p - z*d */
1427: MatMult(jac->H,ilinkA->y,Hu); /* ||u||_H = u'*H*u */
1428: VecDot(Hu,ilinkA->y,&nrmz2);
1430: /* Compute Lower Bound estimate */
1431: if (iterGKB > jac->gkbdelay) {
1432: lowbnd = 0.0;
1433: for (j=0; j<jac->gkbdelay; j++) {
1434: lowbnd += PetscAbsScalar(vecz[j]*vecz[j]);
1435: }
1436: lowbnd = PetscSqrtReal(lowbnd/PetscAbsScalar(nrmz2));
1437: }
1439: for (j=0; j<jac->gkbdelay-1; j++) {
1440: vecz[jac->gkbdelay-j-1] = vecz[jac->gkbdelay-j-2];
1441: }
1442: if (jac->gkbmonitor) {
1443: PetscViewerASCIIPrintf(jac->gkbviewer,"%3D GKB Lower bound estimate %14.12e\n",iterGKB,lowbnd);
1444: }
1445: }
1447: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1448: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1449: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1450: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1452: return(0);
1453: }
1455: #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \
1456: (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
1457: VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
1458: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \
1459: KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \
1460: KSPCheckSolve(ilink->ksp,pc,ilink->x) || \
1461: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \
1462: VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \
1463: VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE))
1465: static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y)
1466: {
1467: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1468: PetscErrorCode ierr;
1469: PC_FieldSplitLink ilink = jac->head;
1470: PetscInt bs;
1473: if (jac->type == PC_COMPOSITE_ADDITIVE) {
1474: if (jac->defaultsplit) {
1475: VecGetBlockSize(x,&bs);
1476: if (jac->bs > 0 && bs != jac->bs) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Blocksize of x vector %D does not match fieldsplit blocksize %D",bs,jac->bs);
1477: VecGetBlockSize(y,&bs);
1478: if (jac->bs > 0 && bs != jac->bs) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Blocksize of y vector %D does not match fieldsplit blocksize %D",bs,jac->bs);
1479: VecStrideGatherAll(x,jac->x,INSERT_VALUES);
1480: while (ilink) {
1481: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1482: KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);
1483: KSPCheckSolve(ilink->ksp,pc,ilink->y);
1484: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1485: ilink = ilink->next;
1486: }
1487: VecStrideScatterAll(jac->y,y,INSERT_VALUES);
1488: } else {
1489: VecSet(y,0.0);
1490: while (ilink) {
1491: FieldSplitSplitSolveAddTranspose(ilink,x,y);
1492: ilink = ilink->next;
1493: }
1494: }
1495: } else {
1496: if (!jac->w1) {
1497: VecDuplicate(x,&jac->w1);
1498: VecDuplicate(x,&jac->w2);
1499: }
1500: VecSet(y,0.0);
1501: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1502: FieldSplitSplitSolveAddTranspose(ilink,x,y);
1503: while (ilink->next) {
1504: ilink = ilink->next;
1505: MatMultTranspose(pc->mat,y,jac->w1);
1506: VecWAXPY(jac->w2,-1.0,jac->w1,x);
1507: FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1508: }
1509: while (ilink->previous) {
1510: ilink = ilink->previous;
1511: MatMultTranspose(pc->mat,y,jac->w1);
1512: VecWAXPY(jac->w2,-1.0,jac->w1,x);
1513: FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1514: }
1515: } else {
1516: while (ilink->next) { /* get to last entry in linked list */
1517: ilink = ilink->next;
1518: }
1519: FieldSplitSplitSolveAddTranspose(ilink,x,y);
1520: while (ilink->previous) {
1521: ilink = ilink->previous;
1522: MatMultTranspose(pc->mat,y,jac->w1);
1523: VecWAXPY(jac->w2,-1.0,jac->w1,x);
1524: FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1525: }
1526: }
1527: }
1528: return(0);
1529: }
1531: static PetscErrorCode PCReset_FieldSplit(PC pc)
1532: {
1533: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1534: PetscErrorCode ierr;
1535: PC_FieldSplitLink ilink = jac->head,next;
1538: while (ilink) {
1539: KSPDestroy(&ilink->ksp);
1540: VecDestroy(&ilink->x);
1541: VecDestroy(&ilink->y);
1542: VecDestroy(&ilink->z);
1543: VecScatterDestroy(&ilink->sctx);
1544: ISDestroy(&ilink->is);
1545: ISDestroy(&ilink->is_col);
1546: PetscFree(ilink->splitname);
1547: PetscFree(ilink->fields);
1548: PetscFree(ilink->fields_col);
1549: next = ilink->next;
1550: PetscFree(ilink);
1551: ilink = next;
1552: }
1553: jac->head = NULL;
1554: PetscFree2(jac->x,jac->y);
1555: if (jac->mat && jac->mat != jac->pmat) {
1556: MatDestroyMatrices(jac->nsplits,&jac->mat);
1557: } else if (jac->mat) {
1558: jac->mat = NULL;
1559: }
1560: if (jac->pmat) {MatDestroyMatrices(jac->nsplits,&jac->pmat);}
1561: if (jac->Afield) {MatDestroyMatrices(jac->nsplits,&jac->Afield);}
1562: jac->nsplits = 0;
1563: VecDestroy(&jac->w1);
1564: VecDestroy(&jac->w2);
1565: MatDestroy(&jac->schur);
1566: MatDestroy(&jac->schurp);
1567: MatDestroy(&jac->schur_user);
1568: KSPDestroy(&jac->kspschur);
1569: KSPDestroy(&jac->kspupper);
1570: MatDestroy(&jac->B);
1571: MatDestroy(&jac->C);
1572: MatDestroy(&jac->H);
1573: VecDestroy(&jac->u);
1574: VecDestroy(&jac->v);
1575: VecDestroy(&jac->Hu);
1576: VecDestroy(&jac->d);
1577: PetscFree(jac->vecz);
1578: PetscViewerDestroy(&jac->gkbviewer);
1579: jac->isrestrict = PETSC_FALSE;
1580: return(0);
1581: }
1583: static PetscErrorCode PCDestroy_FieldSplit(PC pc)
1584: {
1585: PetscErrorCode ierr;
1588: PCReset_FieldSplit(pc);
1589: PetscFree(pc->data);
1590: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurGetSubKSP_C",NULL);
1591: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL);
1592: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL);
1593: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL);
1594: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL);
1595: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL);
1596: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL);
1597: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL);
1598: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL);
1599: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",NULL);
1600: return(0);
1601: }
1603: static PetscErrorCode PCSetFromOptions_FieldSplit(PetscOptionItems *PetscOptionsObject,PC pc)
1604: {
1605: PetscErrorCode ierr;
1606: PetscInt bs;
1607: PetscBool flg;
1608: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1609: PCCompositeType ctype;
1612: PetscOptionsHead(PetscOptionsObject,"FieldSplit options");
1613: PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL);
1614: PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);
1615: if (flg) {
1616: PCFieldSplitSetBlockSize(pc,bs);
1617: }
1618: jac->diag_use_amat = pc->useAmat;
1619: 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);
1620: jac->offdiag_use_amat = pc->useAmat;
1621: 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);
1622: PetscOptionsBool("-pc_fieldsplit_detect_saddle_point","Form 2-way split by detecting zero diagonal entries", "PCFieldSplitSetDetectSaddlePoint",jac->detect,&jac->detect,NULL);
1623: PCFieldSplitSetDetectSaddlePoint(pc,jac->detect); /* Sets split type and Schur PC type */
1624: PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);
1625: if (flg) {
1626: PCFieldSplitSetType(pc,ctype);
1627: }
1628: /* Only setup fields once */
1629: if ((jac->bs > 0) && (jac->nsplits == 0)) {
1630: /* only allow user to set fields from command line if bs is already known.
1631: otherwise user can set them in PCFieldSplitSetDefaults() */
1632: PCFieldSplitSetRuntimeSplits_Private(pc);
1633: if (jac->splitdefined) {PetscInfo(pc,"Splits defined using the options database\n");}
1634: }
1635: if (jac->type == PC_COMPOSITE_SCHUR) {
1636: PetscOptionsGetEnum(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);
1637: if (flg) {PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");}
1638: 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);
1639: PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSetSchurPre",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL);
1640: PetscOptionsScalar("-pc_fieldsplit_schur_scale","Scale Schur complement","PCFieldSplitSetSchurScale",jac->schurscale,&jac->schurscale,NULL);
1641: } else if (jac->type == PC_COMPOSITE_GKB) {
1642: PetscOptionsReal("-pc_fieldsplit_gkb_tol","The tolerance for the lower bound stopping criterion","PCFieldSplitGKBTol",jac->gkbtol,&jac->gkbtol,NULL);
1643: PetscOptionsInt("-pc_fieldsplit_gkb_delay","The delay value for lower bound criterion","PCFieldSplitGKBDelay",jac->gkbdelay,&jac->gkbdelay,NULL);
1644: PetscOptionsReal("-pc_fieldsplit_gkb_nu","Parameter in augmented Lagrangian approach","PCFieldSplitGKBNu",jac->gkbnu,&jac->gkbnu,NULL);
1645: if (jac->gkbnu < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nu cannot be less than 0: value %f",jac->gkbnu);
1646: PetscOptionsInt("-pc_fieldsplit_gkb_maxit","Maximum allowed number of iterations","PCFieldSplitGKBMaxit",jac->gkbmaxit,&jac->gkbmaxit,NULL);
1647: PetscOptionsBool("-pc_fieldsplit_gkb_monitor","Prints number of GKB iterations and error","PCFieldSplitGKB",jac->gkbmonitor,&jac->gkbmonitor,NULL);
1648: }
1649: PetscOptionsTail();
1650: return(0);
1651: }
1653: /*------------------------------------------------------------------------------------*/
1655: static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1656: {
1657: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1658: PetscErrorCode ierr;
1659: PC_FieldSplitLink ilink,next = jac->head;
1660: char prefix[128];
1661: PetscInt i;
1664: if (jac->splitdefined) {
1665: PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1666: return(0);
1667: }
1668: for (i=0; i<n; i++) {
1669: if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs);
1670: if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]);
1671: }
1672: PetscNew(&ilink);
1673: if (splitname) {
1674: PetscStrallocpy(splitname,&ilink->splitname);
1675: } else {
1676: PetscMalloc1(3,&ilink->splitname);
1677: PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);
1678: }
1679: 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 */
1680: PetscMalloc1(n,&ilink->fields);
1681: PetscArraycpy(ilink->fields,fields,n);
1682: PetscMalloc1(n,&ilink->fields_col);
1683: PetscArraycpy(ilink->fields_col,fields_col,n);
1685: ilink->nfields = n;
1686: ilink->next = NULL;
1687: KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1688: KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);
1689: PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1690: KSPSetType(ilink->ksp,KSPPREONLY);
1691: PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);
1693: PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);
1694: KSPSetOptionsPrefix(ilink->ksp,prefix);
1696: if (!next) {
1697: jac->head = ilink;
1698: ilink->previous = NULL;
1699: } else {
1700: while (next->next) {
1701: next = next->next;
1702: }
1703: next->next = ilink;
1704: ilink->previous = next;
1705: }
1706: jac->nsplits++;
1707: return(0);
1708: }
1710: static PetscErrorCode PCFieldSplitSchurGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp)
1711: {
1712: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1716: *subksp = NULL;
1717: if (n) *n = 0;
1718: if (jac->type == PC_COMPOSITE_SCHUR) {
1719: PetscInt nn;
1721: if (!jac->schur) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitSchurGetSubKSP()");
1722: if (jac->nsplits != 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Unexpected number of splits %D != 2",jac->nsplits);
1723: nn = jac->nsplits + (jac->kspupper != jac->head->ksp ? 1 : 0);
1724: PetscMalloc1(nn,subksp);
1725: (*subksp)[0] = jac->head->ksp;
1726: (*subksp)[1] = jac->kspschur;
1727: if (jac->kspupper != jac->head->ksp) (*subksp)[2] = jac->kspupper;
1728: if (n) *n = nn;
1729: }
1730: return(0);
1731: }
1733: static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp)
1734: {
1735: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1739: if (!jac->schur) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitGetSubKSP()");
1740: PetscMalloc1(jac->nsplits,subksp);
1741: MatSchurComplementGetKSP(jac->schur,*subksp);
1743: (*subksp)[1] = jac->kspschur;
1744: if (n) *n = jac->nsplits;
1745: return(0);
1746: }
1748: static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp)
1749: {
1750: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1751: PetscErrorCode ierr;
1752: PetscInt cnt = 0;
1753: PC_FieldSplitLink ilink = jac->head;
1756: PetscMalloc1(jac->nsplits,subksp);
1757: while (ilink) {
1758: (*subksp)[cnt++] = ilink->ksp;
1759: ilink = ilink->next;
1760: }
1761: if (cnt != jac->nsplits) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupt PCFIELDSPLIT object: number of splits in linked list %D does not match number in object %D",cnt,jac->nsplits);
1762: if (n) *n = jac->nsplits;
1763: return(0);
1764: }
1766: /*@C
1767: PCFieldSplitRestrictIS - Restricts the fieldsplit ISs to be within a given IS.
1769: Input Parameters:
1770: + pc - the preconditioner context
1771: - is - the index set that defines the indices to which the fieldsplit is to be restricted
1773: Level: advanced
1775: @*/
1776: PetscErrorCode PCFieldSplitRestrictIS(PC pc,IS isy)
1777: {
1783: PetscTryMethod(pc,"PCFieldSplitRestrictIS_C",(PC,IS),(pc,isy));
1784: return(0);
1785: }
1787: static PetscErrorCode PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy)
1788: {
1789: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1790: PetscErrorCode ierr;
1791: PC_FieldSplitLink ilink = jac->head, next;
1792: PetscInt localsize,size,sizez,i;
1793: const PetscInt *ind, *indz;
1794: PetscInt *indc, *indcz;
1795: PetscBool flg;
1798: ISGetLocalSize(isy,&localsize);
1799: MPI_Scan(&localsize,&size,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)isy));
1800: size -= localsize;
1801: while (ilink) {
1802: IS isrl,isr;
1803: PC subpc;
1804: ISEmbed(ilink->is, isy, PETSC_TRUE, &isrl);
1805: ISGetLocalSize(isrl,&localsize);
1806: PetscMalloc1(localsize,&indc);
1807: ISGetIndices(isrl,&ind);
1808: PetscArraycpy(indc,ind,localsize);
1809: ISRestoreIndices(isrl,&ind);
1810: ISDestroy(&isrl);
1811: for (i=0; i<localsize; i++) *(indc+i) += size;
1812: ISCreateGeneral(PetscObjectComm((PetscObject)isy),localsize,indc,PETSC_OWN_POINTER,&isr);
1813: PetscObjectReference((PetscObject)isr);
1814: ISDestroy(&ilink->is);
1815: ilink->is = isr;
1816: PetscObjectReference((PetscObject)isr);
1817: ISDestroy(&ilink->is_col);
1818: ilink->is_col = isr;
1819: ISDestroy(&isr);
1820: KSPGetPC(ilink->ksp, &subpc);
1821: PetscObjectTypeCompare((PetscObject)subpc,PCFIELDSPLIT,&flg);
1822: if (flg) {
1823: IS iszl,isz;
1824: MPI_Comm comm;
1825: ISGetLocalSize(ilink->is,&localsize);
1826: comm = PetscObjectComm((PetscObject)ilink->is);
1827: ISEmbed(isy, ilink->is, PETSC_TRUE, &iszl);
1828: MPI_Scan(&localsize,&sizez,1,MPIU_INT,MPI_SUM,comm);
1829: sizez -= localsize;
1830: ISGetLocalSize(iszl,&localsize);
1831: PetscMalloc1(localsize,&indcz);
1832: ISGetIndices(iszl,&indz);
1833: PetscArraycpy(indcz,indz,localsize);
1834: ISRestoreIndices(iszl,&indz);
1835: ISDestroy(&iszl);
1836: for (i=0; i<localsize; i++) *(indcz+i) += sizez;
1837: ISCreateGeneral(comm,localsize,indcz,PETSC_OWN_POINTER,&isz);
1838: PCFieldSplitRestrictIS(subpc,isz);
1839: ISDestroy(&isz);
1840: }
1841: next = ilink->next;
1842: ilink = next;
1843: }
1844: jac->isrestrict = PETSC_TRUE;
1845: return(0);
1846: }
1848: static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is)
1849: {
1850: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1851: PetscErrorCode ierr;
1852: PC_FieldSplitLink ilink, next = jac->head;
1853: char prefix[128];
1856: if (jac->splitdefined) {
1857: PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1858: return(0);
1859: }
1860: PetscNew(&ilink);
1861: if (splitname) {
1862: PetscStrallocpy(splitname,&ilink->splitname);
1863: } else {
1864: PetscMalloc1(8,&ilink->splitname);
1865: PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);
1866: }
1867: 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 */
1868: PetscObjectReference((PetscObject)is);
1869: ISDestroy(&ilink->is);
1870: ilink->is = is;
1871: PetscObjectReference((PetscObject)is);
1872: ISDestroy(&ilink->is_col);
1873: ilink->is_col = is;
1874: ilink->next = NULL;
1875: KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1876: KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);
1877: PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1878: KSPSetType(ilink->ksp,KSPPREONLY);
1879: PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);
1881: PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);
1882: KSPSetOptionsPrefix(ilink->ksp,prefix);
1884: if (!next) {
1885: jac->head = ilink;
1886: ilink->previous = NULL;
1887: } else {
1888: while (next->next) {
1889: next = next->next;
1890: }
1891: next->next = ilink;
1892: ilink->previous = next;
1893: }
1894: jac->nsplits++;
1895: return(0);
1896: }
1898: /*@C
1899: PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner
1901: Logically Collective on PC
1903: Input Parameters:
1904: + pc - the preconditioner context
1905: . splitname - name of this split, if NULL the number of the split is used
1906: . n - the number of fields in this split
1907: - fields - the fields in this split
1909: Level: intermediate
1911: Notes:
1912: Use PCFieldSplitSetIS() to set a completely general set of indices as a field.
1914: The PCFieldSplitSetFields() is for defining fields as strided blocks. For example, if the block
1915: 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
1916: 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x....
1917: where the numbered entries indicate what is in the field.
1919: This function is called once per split (it creates a new split each time). Solve options
1920: for this split will be available under the prefix -fieldsplit_SPLITNAME_.
1922: Developer Note: This routine does not actually create the IS representing the split, that is delayed
1923: until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be
1924: available when this routine is called.
1926: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS()
1928: @*/
1929: PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1930: {
1936: if (n < 1) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname);
1938: PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt*,const PetscInt*),(pc,splitname,n,fields,fields_col));
1939: return(0);
1940: }
1942: /*@
1943: PCFieldSplitSetDiagUseAmat - set flag indicating whether to extract diagonal blocks from Amat (rather than Pmat)
1945: Logically Collective on PC
1947: Input Parameters:
1948: + pc - the preconditioner object
1949: - flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from
1951: Options Database:
1952: . -pc_fieldsplit_diag_use_amat
1954: Level: intermediate
1956: .seealso: PCFieldSplitGetDiagUseAmat(), PCFieldSplitSetOffDiagUseAmat(), PCFIELDSPLIT
1958: @*/
1959: PetscErrorCode PCFieldSplitSetDiagUseAmat(PC pc,PetscBool flg)
1960: {
1961: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1962: PetscBool isfs;
1967: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1968: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1969: jac->diag_use_amat = flg;
1970: return(0);
1971: }
1973: /*@
1974: PCFieldSplitGetDiagUseAmat - get the flag indicating whether to extract diagonal blocks from Amat (rather than Pmat)
1976: Logically Collective on PC
1978: Input Parameters:
1979: . pc - the preconditioner object
1981: Output Parameters:
1982: . flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from
1984: Level: intermediate
1986: .seealso: PCFieldSplitSetDiagUseAmat(), PCFieldSplitGetOffDiagUseAmat(), PCFIELDSPLIT
1988: @*/
1989: PetscErrorCode PCFieldSplitGetDiagUseAmat(PC pc,PetscBool *flg)
1990: {
1991: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1992: PetscBool isfs;
1998: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1999: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
2000: *flg = jac->diag_use_amat;
2001: return(0);
2002: }
2004: /*@
2005: PCFieldSplitSetOffDiagUseAmat - set flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat)
2007: Logically Collective on PC
2009: Input Parameters:
2010: + pc - the preconditioner object
2011: - flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from
2013: Options Database:
2014: . -pc_fieldsplit_off_diag_use_amat
2016: Level: intermediate
2018: .seealso: PCFieldSplitGetOffDiagUseAmat(), PCFieldSplitSetDiagUseAmat(), PCFIELDSPLIT
2020: @*/
2021: PetscErrorCode PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg)
2022: {
2023: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2024: PetscBool isfs;
2029: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2030: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
2031: jac->offdiag_use_amat = flg;
2032: return(0);
2033: }
2035: /*@
2036: PCFieldSplitGetOffDiagUseAmat - get the flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat)
2038: Logically Collective on PC
2040: Input Parameters:
2041: . pc - the preconditioner object
2043: Output Parameters:
2044: . flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from
2046: Level: intermediate
2048: .seealso: PCFieldSplitSetOffDiagUseAmat(), PCFieldSplitGetDiagUseAmat(), PCFIELDSPLIT
2050: @*/
2051: PetscErrorCode PCFieldSplitGetOffDiagUseAmat(PC pc,PetscBool *flg)
2052: {
2053: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2054: PetscBool isfs;
2060: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2061: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
2062: *flg = jac->offdiag_use_amat;
2063: return(0);
2064: }
2066: /*@C
2067: PCFieldSplitSetIS - Sets the exact elements for field
2069: Logically Collective on PC
2071: Input Parameters:
2072: + pc - the preconditioner context
2073: . splitname - name of this split, if NULL the number of the split is used
2074: - is - the index set that defines the vector elements in this field
2076: Notes:
2077: Use PCFieldSplitSetFields(), for fields defined by strided types.
2079: This function is called once per split (it creates a new split each time). Solve options
2080: for this split will be available under the prefix -fieldsplit_SPLITNAME_.
2082: Level: intermediate
2084: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize()
2086: @*/
2087: PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is)
2088: {
2095: PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));
2096: return(0);
2097: }
2099: /*@C
2100: PCFieldSplitGetIS - Retrieves the elements for a field as an IS
2102: Logically Collective on PC
2104: Input Parameters:
2105: + pc - the preconditioner context
2106: - splitname - name of this split
2108: Output Parameter:
2109: - is - the index set that defines the vector elements in this field, or NULL if the field is not found
2111: Level: intermediate
2113: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS()
2115: @*/
2116: PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is)
2117: {
2124: {
2125: PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;
2126: PC_FieldSplitLink ilink = jac->head;
2127: PetscBool found;
2129: *is = NULL;
2130: while (ilink) {
2131: PetscStrcmp(ilink->splitname, splitname, &found);
2132: if (found) {
2133: *is = ilink->is;
2134: break;
2135: }
2136: ilink = ilink->next;
2137: }
2138: }
2139: return(0);
2140: }
2142: /*@C
2143: PCFieldSplitGetISByIndex - Retrieves the elements for a given index field as an IS
2145: Logically Collective on PC
2147: Input Parameters:
2148: + pc - the preconditioner context
2149: - index - index of this split
2151: Output Parameter:
2152: - is - the index set that defines the vector elements in this field
2154: Level: intermediate
2156: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitGetIS(), PCFieldSplitSetIS()
2158: @*/
2159: PetscErrorCode PCFieldSplitGetISByIndex(PC pc,PetscInt index,IS *is)
2160: {
2164: if (index < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",index);
2167: {
2168: PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;
2169: PC_FieldSplitLink ilink = jac->head;
2170: PetscInt i = 0;
2171: if (index >= jac->nsplits) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",index,jac->nsplits);
2173: while (i < index) {
2174: ilink = ilink->next;
2175: ++i;
2176: }
2177: PCFieldSplitGetIS(pc,ilink->splitname,is);
2178: }
2179: return(0);
2180: }
2182: /*@
2183: PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the
2184: fieldsplit preconditioner. If not set the matrix block size is used.
2186: Logically Collective on PC
2188: Input Parameters:
2189: + pc - the preconditioner context
2190: - bs - the block size
2192: Level: intermediate
2194: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields()
2196: @*/
2197: PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs)
2198: {
2204: PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));
2205: return(0);
2206: }
2208: /*@C
2209: PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits
2211: Collective on KSP
2213: Input Parameter:
2214: . pc - the preconditioner context
2216: Output Parameters:
2217: + n - the number of splits
2218: - subksp - the array of KSP contexts
2220: Note:
2221: After PCFieldSplitGetSubKSP() the array of KSPs is to be freed by the user with PetscFree()
2222: (not the KSP just the array that contains them).
2224: You must call PCSetUp() before calling PCFieldSplitGetSubKSP().
2226: If the fieldsplit is of type PC_COMPOSITE_SCHUR, it returns the KSP object used inside the
2227: Schur complement and the KSP object used to iterate over the Schur complement.
2228: To access all the KSP objects used in PC_COMPOSITE_SCHUR, use PCFieldSplitSchurGetSubKSP().
2230: If the fieldsplit is of type PC_COMPOSITE_GKB, it returns the KSP object used to solve the
2231: inner linear system defined by the matrix H in each loop.
2233: Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs.
2234: You can call PCFieldSplitGetSubKSP(pc,n,PETSC_NULL_KSP,ierr) to determine how large the
2235: KSP array must be.
2237: Level: advanced
2239: .seealso: PCFIELDSPLIT
2240: @*/
2241: PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
2242: {
2248: PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));
2249: return(0);
2250: }
2252: /*@C
2253: PCFieldSplitSchurGetSubKSP - Gets the KSP contexts used inside the Schur complement based PCFIELDSPLIT
2255: Collective on KSP
2257: Input Parameter:
2258: . pc - the preconditioner context
2260: Output Parameters:
2261: + n - the number of splits
2262: - subksp - the array of KSP contexts
2264: Note:
2265: After PCFieldSplitSchurGetSubKSP() the array of KSPs is to be freed by the user with PetscFree()
2266: (not the KSP just the array that contains them).
2268: You must call PCSetUp() before calling PCFieldSplitSchurGetSubKSP().
2270: If the fieldsplit type is of type PC_COMPOSITE_SCHUR, it returns (in order)
2271: - the KSP used for the (1,1) block
2272: - the KSP used for the Schur complement (not the one used for the interior Schur solver)
2273: - the KSP used for the (1,1) block in the upper triangular factor (if different from that of the (1,1) block).
2275: It returns a null array if the fieldsplit is not of type PC_COMPOSITE_SCHUR; in this case, you should use PCFieldSplitGetSubKSP().
2277: Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs.
2278: You can call PCFieldSplitSchurGetSubKSP(pc,n,PETSC_NULL_KSP,ierr) to determine how large the
2279: KSP array must be.
2281: Level: advanced
2283: .seealso: PCFIELDSPLIT
2284: @*/
2285: PetscErrorCode PCFieldSplitSchurGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
2286: {
2292: PetscUseMethod(pc,"PCFieldSplitSchurGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));
2293: return(0);
2294: }
2296: /*@
2297: PCFieldSplitSetSchurPre - Indicates from what operator the preconditioner is constructucted for the Schur complement.
2298: The default is the A11 matrix.
2300: Collective on PC
2302: Input Parameters:
2303: + pc - the preconditioner context
2304: . 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
2305: PC_FIELDSPLIT_SCHUR_PRE_SELFP, and PC_FIELDSPLIT_SCHUR_PRE_FULL
2306: - userpre - matrix to use for preconditioning, or NULL
2308: Options Database:
2309: + -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11. See notes for meaning of various arguments
2310: - -fieldsplit_1_pc_type <pctype> - the preconditioner algorithm that is used to construct the preconditioner from the operator
2312: Notes:
2313: $ If ptype is
2314: $ a11 - the preconditioner for the Schur complement is generated from the block diagonal part of the preconditioner
2315: $ matrix associated with the Schur complement (i.e. A11), not the Schur complement matrix
2316: $ self - the preconditioner for the Schur complement is generated from the symbolic representation of the Schur complement matrix:
2317: $ The only preconditioner that currently works with this symbolic respresentation matrix object is the PCLSC
2318: $ preconditioner
2319: $ user - the preconditioner for the Schur complement is generated from the user provided matrix (pre argument
2320: $ to this function).
2321: $ selfp - the preconditioning for the Schur complement is generated from an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01
2322: $ This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be
2323: $ lumped before extracting the diagonal using the additional option -fieldsplit_1_mat_schur_complement_ainv_type lump
2324: $ full - the preconditioner for the Schur complement is generated from the exact Schur complement matrix representation computed internally by PCFIELDSPLIT (this is expensive)
2325: $ useful mostly as a test that the Schur complement approach can work for your problem
2327: When solving a saddle point problem, where the A11 block is identically zero, using a11 as the ptype only makes sense
2328: with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and
2329: -fieldsplit_1_pc_type lsc which uses the least squares commutator to compute a preconditioner for the Schur complement.
2331: Level: intermediate
2333: .seealso: PCFieldSplitGetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType,
2334: MatSchurComplementSetAinvType(), PCLSC
2336: @*/
2337: PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
2338: {
2343: PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));
2344: return(0);
2345: }
2347: PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) {return PCFieldSplitSetSchurPre(pc,ptype,pre);} /* Deprecated name */
2349: /*@
2350: PCFieldSplitGetSchurPre - For Schur complement fieldsplit, determine how the Schur complement will be
2351: preconditioned. See PCFieldSplitSetSchurPre() for details.
2353: Logically Collective on PC
2355: Input Parameter:
2356: . pc - the preconditioner context
2358: Output Parameters:
2359: + ptype - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11, PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_PRE_USER
2360: - userpre - matrix to use for preconditioning (with PC_FIELDSPLIT_PRE_USER), or NULL
2362: Level: intermediate
2364: .seealso: PCFieldSplitSetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC
2366: @*/
2367: PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
2368: {
2373: PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre));
2374: return(0);
2375: }
2377: /*@
2378: PCFieldSplitSchurGetS - extract the MatSchurComplement object used by this PC in case it needs to be configured separately
2380: Not collective
2382: Input Parameter:
2383: . pc - the preconditioner context
2385: Output Parameter:
2386: . S - the Schur complement matrix
2388: Notes:
2389: This matrix should not be destroyed using MatDestroy(); rather, use PCFieldSplitSchurRestoreS().
2391: Level: advanced
2393: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurRestoreS()
2395: @*/
2396: PetscErrorCode PCFieldSplitSchurGetS(PC pc,Mat *S)
2397: {
2399: const char* t;
2400: PetscBool isfs;
2401: PC_FieldSplit *jac;
2405: PetscObjectGetType((PetscObject)pc,&t);
2406: PetscStrcmp(t,PCFIELDSPLIT,&isfs);
2407: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
2408: jac = (PC_FieldSplit*)pc->data;
2409: if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
2410: if (S) *S = jac->schur;
2411: return(0);
2412: }
2414: /*@
2415: PCFieldSplitSchurRestoreS - restores the MatSchurComplement object used by this PC
2417: Not collective
2419: Input Parameters:
2420: + pc - the preconditioner context
2421: - S - the Schur complement matrix
2423: Level: advanced
2425: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurGetS()
2427: @*/
2428: PetscErrorCode PCFieldSplitSchurRestoreS(PC pc,Mat *S)
2429: {
2431: const char* t;
2432: PetscBool isfs;
2433: PC_FieldSplit *jac;
2437: PetscObjectGetType((PetscObject)pc,&t);
2438: PetscStrcmp(t,PCFIELDSPLIT,&isfs);
2439: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
2440: jac = (PC_FieldSplit*)pc->data;
2441: if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
2442: if (!S || *S != jac->schur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten");
2443: return(0);
2444: }
2446: static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
2447: {
2448: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2452: jac->schurpre = ptype;
2453: if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) {
2454: MatDestroy(&jac->schur_user);
2455: jac->schur_user = pre;
2456: PetscObjectReference((PetscObject)jac->schur_user);
2457: }
2458: return(0);
2459: }
2461: static PetscErrorCode PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
2462: {
2463: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2466: *ptype = jac->schurpre;
2467: *pre = jac->schur_user;
2468: return(0);
2469: }
2471: /*@
2472: PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain in the preconditioner
2474: Collective on PC
2476: Input Parameters:
2477: + pc - the preconditioner context
2478: - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default
2480: Options Database:
2481: . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full
2483: Level: intermediate
2485: Notes:
2486: The FULL factorization is
2488: $ (A B) = (1 0) (A 0) (1 Ainv*B) = L D U
2489: $ (C E) (C*Ainv 1) (0 S) (0 1)
2491: 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,
2492: 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().
2494: $ If A and S are solved exactly
2495: $ *) FULL factorization is a direct solver.
2496: $ *) 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.
2497: $ *) 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.
2499: If the iteration count is very low, consider using KSPFGMRES or KSPGCR which can use one less preconditioner
2500: application in this case. Note that the preconditioned operator may be highly non-normal, so such fast convergence may not be observed in practice.
2502: For symmetric problems in which A is positive definite and S is negative definite, DIAG can be used with KSPMINRES.
2504: 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).
2506: References:
2507: + 1. - Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000).
2508: - 2. - Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001).
2510: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCFieldSplitSetSchurScale()
2511: @*/
2512: PetscErrorCode PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype)
2513: {
2518: PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));
2519: return(0);
2520: }
2522: static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype)
2523: {
2524: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2527: jac->schurfactorization = ftype;
2528: return(0);
2529: }
2531: /*@
2532: PCFieldSplitSetSchurScale - Controls the sign flip of S for PC_FIELDSPLIT_SCHUR_FACT_DIAG.
2534: Collective on PC
2536: Input Parameters:
2537: + pc - the preconditioner context
2538: - scale - scaling factor for the Schur complement
2540: Options Database:
2541: . -pc_fieldsplit_schur_scale - default is -1.0
2543: Level: intermediate
2545: .seealso: PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurFactType, PCFieldSplitSetSchurScale()
2546: @*/
2547: PetscErrorCode PCFieldSplitSetSchurScale(PC pc,PetscScalar scale)
2548: {
2554: PetscTryMethod(pc,"PCFieldSplitSetSchurScale_C",(PC,PetscScalar),(pc,scale));
2555: return(0);
2556: }
2558: static PetscErrorCode PCFieldSplitSetSchurScale_FieldSplit(PC pc,PetscScalar scale)
2559: {
2560: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2563: jac->schurscale = scale;
2564: return(0);
2565: }
2567: /*@C
2568: PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement
2570: Collective on KSP
2572: Input Parameter:
2573: . pc - the preconditioner context
2575: Output Parameters:
2576: + A00 - the (0,0) block
2577: . A01 - the (0,1) block
2578: . A10 - the (1,0) block
2579: - A11 - the (1,1) block
2581: Level: advanced
2583: .seealso: PCFIELDSPLIT
2584: @*/
2585: PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11)
2586: {
2587: PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;
2591: if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach.");
2592: if (A00) *A00 = jac->pmat[0];
2593: if (A01) *A01 = jac->B;
2594: if (A10) *A10 = jac->C;
2595: if (A11) *A11 = jac->pmat[1];
2596: return(0);
2597: }
2599: /*@
2600: PCFieldSplitSetGKBTol - Sets the solver tolerance for the generalized Golub-Kahan bidiagonalization preconditioner.
2602: Collective on PC
2604: Notes:
2605: The generalized GKB algorithm uses a lower bound estimate of the error in energy norm as stopping criterion.
2606: It stops once the lower bound estimate undershoots the required solver tolerance. Although the actual error might be bigger than
2607: this estimate, the stopping criterion is satisfactory in practical cases [A13].
2609: [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013.
2611: Input Parameters:
2612: + pc - the preconditioner context
2613: - tolerance - the solver tolerance
2615: Options Database:
2616: . -pc_fieldsplit_gkb_tol - default is 1e-5
2618: Level: intermediate
2620: .seealso: PCFIELDSPLIT, PCFieldSplitSetGKBDelay(), PCFieldSplitSetGKBNu(), PCFieldSplitSetGKBMaxit()
2621: @*/
2622: PetscErrorCode PCFieldSplitSetGKBTol(PC pc,PetscReal tolerance)
2623: {
2629: PetscTryMethod(pc,"PCFieldSplitSetGKBTol_C",(PC,PetscReal),(pc,tolerance));
2630: return(0);
2631: }
2633: static PetscErrorCode PCFieldSplitSetGKBTol_FieldSplit(PC pc,PetscReal tolerance)
2634: {
2635: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2638: jac->gkbtol = tolerance;
2639: return(0);
2640: }
2642: /*@
2643: PCFieldSplitSetGKBMaxit - Sets the maximum number of iterations for the generalized Golub-Kahan bidiagonalization
2644: preconditioner.
2646: Collective on PC
2648: Input Parameters:
2649: + pc - the preconditioner context
2650: - maxit - the maximum number of iterations
2652: Options Database:
2653: . -pc_fieldsplit_gkb_maxit - default is 100
2655: Level: intermediate
2657: .seealso: PCFIELDSPLIT, PCFieldSplitSetGKBDelay(), PCFieldSplitSetGKBTol(), PCFieldSplitSetGKBNu()
2658: @*/
2659: PetscErrorCode PCFieldSplitSetGKBMaxit(PC pc,PetscInt maxit)
2660: {
2666: PetscTryMethod(pc,"PCFieldSplitSetGKBMaxit_C",(PC,PetscInt),(pc,maxit));
2667: return(0);
2668: }
2670: static PetscErrorCode PCFieldSplitSetGKBMaxit_FieldSplit(PC pc,PetscInt maxit)
2671: {
2672: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2675: jac->gkbmaxit = maxit;
2676: return(0);
2677: }
2679: /*@
2680: PCFieldSplitSetGKBDelay - Sets the delay in the lower bound error estimate in the generalized Golub-Kahan bidiagonalization
2681: preconditioner.
2683: Collective on PC
2685: Notes:
2686: 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
2687: is expressed as a truncated sum. The error at iteration k can only be measured at iteration (k + delay), and thus the algorithm needs
2688: 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
2690: [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013.
2692: Input Parameters:
2693: + pc - the preconditioner context
2694: - delay - the delay window in the lower bound estimate
2696: Options Database:
2697: . -pc_fieldsplit_gkb_delay - default is 5
2699: Level: intermediate
2701: .seealso: PCFIELDSPLIT, PCFieldSplitSetGKBNu(), PCFieldSplitSetGKBTol(), PCFieldSplitSetGKBMaxit()
2702: @*/
2703: PetscErrorCode PCFieldSplitSetGKBDelay(PC pc,PetscInt delay)
2704: {
2710: PetscTryMethod(pc,"PCFieldSplitSetGKBDelay_C",(PC,PetscInt),(pc,delay));
2711: return(0);
2712: }
2714: static PetscErrorCode PCFieldSplitSetGKBDelay_FieldSplit(PC pc,PetscInt delay)
2715: {
2716: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2719: jac->gkbdelay = delay;
2720: return(0);
2721: }
2723: /*@
2724: 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.
2726: Collective on PC
2728: Notes:
2729: 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,
2730: 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
2731: necessary to find a good balance in between the convergence of the inner and outer loop.
2733: 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.
2735: [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013.
2737: Input Parameters:
2738: + pc - the preconditioner context
2739: - nu - the shift parameter
2741: Options Database:
2742: . -pc_fieldsplit_gkb_nu - default is 1
2744: Level: intermediate
2746: .seealso: PCFIELDSPLIT, PCFieldSplitSetGKBDelay(), PCFieldSplitSetGKBTol(), PCFieldSplitSetGKBMaxit()
2747: @*/
2748: PetscErrorCode PCFieldSplitSetGKBNu(PC pc,PetscReal nu)
2749: {
2755: PetscTryMethod(pc,"PCFieldSplitSetGKBNu_C",(PC,PetscReal),(pc,nu));
2756: return(0);
2757: }
2759: static PetscErrorCode PCFieldSplitSetGKBNu_FieldSplit(PC pc,PetscReal nu)
2760: {
2761: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2764: jac->gkbnu = nu;
2765: return(0);
2766: }
2768: static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type)
2769: {
2770: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2774: jac->type = type;
2776: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",0);
2777: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",0);
2778: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",0);
2779: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);
2780: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",0);
2781: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBTol_C",0);
2782: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBMaxit_C",0);
2783: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBNu_C",0);
2784: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBDelay_C",0);
2786: if (type == PC_COMPOSITE_SCHUR) {
2787: pc->ops->apply = PCApply_FieldSplit_Schur;
2788: pc->ops->view = PCView_FieldSplit_Schur;
2790: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);
2791: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit);
2792: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit);
2793: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);
2794: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",PCFieldSplitSetSchurScale_FieldSplit);
2795: } else if (type == PC_COMPOSITE_GKB) {
2796: pc->ops->apply = PCApply_FieldSplit_GKB;
2797: pc->ops->view = PCView_FieldSplit_GKB;
2799: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2800: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBTol_C",PCFieldSplitSetGKBTol_FieldSplit);
2801: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBMaxit_C",PCFieldSplitSetGKBMaxit_FieldSplit);
2802: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBNu_C",PCFieldSplitSetGKBNu_FieldSplit);
2803: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBDelay_C",PCFieldSplitSetGKBDelay_FieldSplit);
2804: } else {
2805: pc->ops->apply = PCApply_FieldSplit;
2806: pc->ops->view = PCView_FieldSplit;
2808: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2809: }
2810: return(0);
2811: }
2813: static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs)
2814: {
2815: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2818: if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs);
2819: if (jac->bs > 0 && jac->bs != bs) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Cannot change fieldsplit blocksize from %D to %D after it has been set",jac->bs,bs);
2820: jac->bs = bs;
2821: return(0);
2822: }
2824: /*@
2825: PCFieldSplitSetType - Sets the type of fieldsplit preconditioner.
2827: Collective on PC
2829: Input Parameters:
2830: + pc - the preconditioner context
2831: - type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR
2833: Options Database Key:
2834: . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type
2836: Level: Intermediate
2838: .seealso: PCCompositeSetType()
2840: @*/
2841: PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type)
2842: {
2847: PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));
2848: return(0);
2849: }
2851: /*@
2852: PCFieldSplitGetType - Gets the type of fieldsplit preconditioner.
2854: Not collective
2856: Input Parameter:
2857: . pc - the preconditioner context
2859: Output Parameter:
2860: . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR
2862: Level: Intermediate
2864: .seealso: PCCompositeSetType()
2865: @*/
2866: PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type)
2867: {
2868: PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;
2873: *type = jac->type;
2874: return(0);
2875: }
2877: /*@
2878: PCFieldSplitSetDMSplits - Flags whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible.
2880: Logically Collective
2882: Input Parameters:
2883: + pc - the preconditioner context
2884: - flg - boolean indicating whether to use field splits defined by the DM
2886: Options Database Key:
2887: . -pc_fieldsplit_dm_splits
2889: Level: Intermediate
2891: .seealso: PCFieldSplitGetDMSplits()
2893: @*/
2894: PetscErrorCode PCFieldSplitSetDMSplits(PC pc,PetscBool flg)
2895: {
2896: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2897: PetscBool isfs;
2903: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2904: if (isfs) {
2905: jac->dm_splits = flg;
2906: }
2907: return(0);
2908: }
2910: /*@
2911: PCFieldSplitGetDMSplits - Returns flag indicating whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible.
2913: Logically Collective
2915: Input Parameter:
2916: . pc - the preconditioner context
2918: Output Parameter:
2919: . flg - boolean indicating whether to use field splits defined by the DM
2921: Level: Intermediate
2923: .seealso: PCFieldSplitSetDMSplits()
2925: @*/
2926: PetscErrorCode PCFieldSplitGetDMSplits(PC pc,PetscBool* flg)
2927: {
2928: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2929: PetscBool isfs;
2935: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2936: if (isfs) {
2937: if (flg) *flg = jac->dm_splits;
2938: }
2939: return(0);
2940: }
2942: /*@
2943: PCFieldSplitGetDetectSaddlePoint - Returns flag indicating whether PCFieldSplit will attempt to automatically determine fields based on zero diagonal entries.
2945: Logically Collective
2947: Input Parameter:
2948: . pc - the preconditioner context
2950: Output Parameter:
2951: . flg - boolean indicating whether to detect fields or not
2953: Level: Intermediate
2955: .seealso: PCFIELDSPLIT, PCFieldSplitSetDetectSaddlePoint()
2957: @*/
2958: PetscErrorCode PCFieldSplitGetDetectSaddlePoint(PC pc,PetscBool *flg)
2959: {
2960: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2963: *flg = jac->detect;
2964: return(0);
2965: }
2967: /*@
2968: PCFieldSplitSetDetectSaddlePoint - Sets flag indicating whether PCFieldSplit will attempt to automatically determine fields based on zero diagonal entries.
2970: Logically Collective
2972: Notes:
2973: Also sets the split type to PC_COMPOSITE_SCHUR (see PCFieldSplitSetType()) and the Schur preconditioner type to PC_FIELDSPLIT_SCHUR_PRE_SELF (see PCFieldSplitSetSchurPre()).
2975: Input Parameter:
2976: . pc - the preconditioner context
2978: Output Parameter:
2979: . flg - boolean indicating whether to detect fields or not
2981: Options Database Key:
2982: . -pc_fieldsplit_detect_saddle_point
2984: Level: Intermediate
2986: .seealso: PCFIELDSPLIT, PCFieldSplitSetDetectSaddlePoint(), PCFieldSplitSetType(), PCFieldSplitSetSchurPre()
2988: @*/
2989: PetscErrorCode PCFieldSplitSetDetectSaddlePoint(PC pc,PetscBool flg)
2990: {
2991: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2995: jac->detect = flg;
2996: if (jac->detect) {
2997: PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);
2998: PCFieldSplitSetSchurPre(pc,PC_FIELDSPLIT_SCHUR_PRE_SELF,NULL);
2999: }
3000: return(0);
3001: }
3003: /* -------------------------------------------------------------------------------------*/
3004: /*MC
3005: PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual
3006: fields or groups of fields. See the users manual section "Solving Block Matrices" for more details.
3008: To set options on the solvers for each block append -fieldsplit_ to all the PC
3009: options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1
3011: To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP()
3012: and set the options directly on the resulting KSP object
3014: Level: intermediate
3016: Options Database Keys:
3017: + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split
3018: . -pc_fieldsplit_default - automatically add any fields to additional splits that have not
3019: been supplied explicitly by -pc_fieldsplit_%d_fields
3020: . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields)
3021: . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur,gkb> - type of relaxation or factorization splitting
3022: . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11; see PCFieldSplitSetSchurPre()
3023: . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> - set factorization type when using -pc_fieldsplit_type schur; see PCFieldSplitSetSchurFactType()
3024: - -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero diagonal and uses Schur complement with no preconditioner as the solver
3026: Options prefixes for inner solvers when using the Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ .
3027: For all other solvers they are -fieldsplit_%d_ for the dth field; use -fieldsplit_ for all fields.
3028: The options prefix for the inner solver when using the Golub-Kahan biadiagonalization preconditioner is -fieldsplit_0_
3030: Notes:
3031: Use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS()
3032: to define a field by an arbitrary collection of entries.
3034: If no fields are set the default is used. The fields are defined by entries strided by bs,
3035: beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(),
3036: if this is not called the block size defaults to the blocksize of the second matrix passed
3037: to KSPSetOperators()/PCSetOperators().
3039: $ For the Schur complement preconditioner if J = [ A00 A01]
3040: $ [ A10 A11]
3041: $ the preconditioner using full factorization is
3042: $ [ I -ksp(A00) A01] [ inv(A00) 0 ] [ I 0]
3043: $ [ 0 I ] [ 0 ksp(S)] [ -A10 ksp(A00) I]
3044: where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. S is the Schur complement
3045: $ S = A11 - A10 ksp(A00) A01
3046: 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
3047: in providing the SECOND split or 1 if not give). For PCFieldSplitGetSubKSP() when field number is 0,
3048: it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default
3049: A11 is used to construct a preconditioner for S, use PCFieldSplitSetSchurPre() for all the possible ways to construct the preconditioner for S.
3051: The factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above,
3052: diag gives
3053: $ [ inv(A00) 0 ]
3054: $ [ 0 -ksp(S)]
3055: 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
3056: can be turned off with PCFieldSplitSetSchurScale() or by command line -pc_fieldsplit_schur_scale 1.0. The lower factorization is the inverse of
3057: $ [ A00 0]
3058: $ [ A10 S]
3059: where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of
3060: $ [ A00 A01]
3061: $ [ 0 S ]
3062: where again the inverses of A00 and S are applied using KSPs.
3064: If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS
3065: is used automatically for a second block.
3067: The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1.
3068: Generally it should be used with the AIJ format.
3070: The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see,
3071: for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT
3072: inside a smoother resulting in "Distributive Smoothers".
3074: References:
3075: A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier-Stokes equations
3076: Howard Elman, V.E. Howle, John Shadid, Robert Shuttleworth, Ray Tuminaro, Journal of Computational Physics 227 (2008) 1790--1808
3077: http://chess.cs.umd.edu/~elman/papers/tax.pdf
3079: The Constrained Pressure Preconditioner (CPR) can be implemented using PCCOMPOSITE with PCGALERKIN. CPR first solves an R A P subsystem, updates the
3080: residual on all variables (PCCompositeSetType(pc,PC_COMPOSITE_MULTIPLICATIVE)), and then applies a simple ILU like preconditioner on all the variables.
3082: The generalized Golub-Kahan bidiagonalization preconditioner (gkb) can be applied to symmetric 2x2 block matrices of the shape
3083: $ [ A00 A01]
3084: $ [ A01' 0 ]
3085: 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'.
3086: 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_.
3088: [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013.
3090: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCLSC,
3091: PCFieldSplitGetSubKSP(), PCFieldSplitSchurGetSubKSP(), PCFieldSplitSetFields(),
3092: PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre(), PCFieldSplitSetSchurFactType(),
3093: MatSchurComplementSetAinvType(), PCFieldSplitSetSchurScale(), PCFieldSplitSetDetectSaddlePoint()
3094: M*/
3096: PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc)
3097: {
3099: PC_FieldSplit *jac;
3102: PetscNewLog(pc,&jac);
3104: jac->bs = -1;
3105: jac->nsplits = 0;
3106: jac->type = PC_COMPOSITE_MULTIPLICATIVE;
3107: jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */
3108: jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL;
3109: jac->schurscale = -1.0;
3110: jac->dm_splits = PETSC_TRUE;
3111: jac->detect = PETSC_FALSE;
3112: jac->gkbtol = 1e-5;
3113: jac->gkbdelay = 5;
3114: jac->gkbnu = 1;
3115: jac->gkbmaxit = 100;
3116: jac->gkbmonitor = PETSC_FALSE;
3118: pc->data = (void*)jac;
3120: pc->ops->apply = PCApply_FieldSplit;
3121: pc->ops->applytranspose = PCApplyTranspose_FieldSplit;
3122: pc->ops->setup = PCSetUp_FieldSplit;
3123: pc->ops->reset = PCReset_FieldSplit;
3124: pc->ops->destroy = PCDestroy_FieldSplit;
3125: pc->ops->setfromoptions = PCSetFromOptions_FieldSplit;
3126: pc->ops->view = PCView_FieldSplit;
3127: pc->ops->applyrichardson = NULL;
3129: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurGetSubKSP_C",PCFieldSplitSchurGetSubKSP_FieldSplit);
3130: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
3131: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);
3132: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);
3133: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);
3134: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);
3135: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",PCFieldSplitRestrictIS_FieldSplit);
3136: return(0);
3137: }