Actual source code: fieldsplit.c
petsc-3.13.6 2020-09-29
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_",0};
6: const char *const PCFieldSplitSchurFactTypes[] = {"DIAG","LOWER","UPPER","FULL","PCFieldSplitSchurFactType","PC_FIELDSPLIT_SCHUR_FACT_",0};
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: * Section 1.5 Writing Application Codes with PETSc-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: }
90: #include <petscdraw.h>
91: static PetscErrorCode PCView_FieldSplit(PC pc,PetscViewer viewer)
92: {
93: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
94: PetscErrorCode ierr;
95: PetscBool iascii,isdraw;
96: PetscInt i,j;
97: PC_FieldSplitLink ilink = jac->head;
100: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
101: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
102: if (iascii) {
103: if (jac->bs > 0) {
104: PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs);
105: } else {
106: PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D\n",PCCompositeTypes[jac->type],jac->nsplits);
107: }
108: if (pc->useAmat) {
109: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");
110: }
111: if (jac->diag_use_amat) {
112: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for diagonal blocks\n");
113: }
114: if (jac->offdiag_use_amat) {
115: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for off-diagonal blocks\n");
116: }
117: PetscViewerASCIIPrintf(viewer," Solver info for each split is in the following KSP objects:\n");
118: for (i=0; i<jac->nsplits; i++) {
119: if (ilink->fields) {
120: PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);
121: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
122: for (j=0; j<ilink->nfields; j++) {
123: if (j > 0) {
124: PetscViewerASCIIPrintf(viewer,",");
125: }
126: PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);
127: }
128: PetscViewerASCIIPrintf(viewer,"\n");
129: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
130: } else {
131: PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);
132: }
133: KSPView(ilink->ksp,viewer);
134: ilink = ilink->next;
135: }
136: }
138: if (isdraw) {
139: PetscDraw draw;
140: PetscReal x,y,w,wd;
142: PetscViewerDrawGetDraw(viewer,0,&draw);
143: PetscDrawGetCurrentPoint(draw,&x,&y);
144: w = 2*PetscMin(1.0 - x,x);
145: wd = w/(jac->nsplits + 1);
146: x = x - wd*(jac->nsplits-1)/2.0;
147: for (i=0; i<jac->nsplits; i++) {
148: PetscDrawPushCurrentPoint(draw,x,y);
149: KSPView(ilink->ksp,viewer);
150: PetscDrawPopCurrentPoint(draw);
151: x += wd;
152: ilink = ilink->next;
153: }
154: }
155: return(0);
156: }
158: static PetscErrorCode PCView_FieldSplit_Schur(PC pc,PetscViewer viewer)
159: {
160: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
161: PetscErrorCode ierr;
162: PetscBool iascii,isdraw;
163: PetscInt i,j;
164: PC_FieldSplitLink ilink = jac->head;
165: MatSchurComplementAinvType atype;
168: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
169: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
170: if (iascii) {
171: if (jac->bs > 0) {
172: PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, blocksize = %D, factorization %s\n",jac->bs,PCFieldSplitSchurFactTypes[jac->schurfactorization]);
173: } else {
174: PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, factorization %s\n",PCFieldSplitSchurFactTypes[jac->schurfactorization]);
175: }
176: if (pc->useAmat) {
177: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");
178: }
179: switch (jac->schurpre) {
180: case PC_FIELDSPLIT_SCHUR_PRE_SELF:
181: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from S itself\n");
182: break;
183: case PC_FIELDSPLIT_SCHUR_PRE_SELFP:
184: MatSchurComplementGetAinvType(jac->schur,&atype);
185: 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;
186: case PC_FIELDSPLIT_SCHUR_PRE_A11:
187: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n");
188: break;
189: case PC_FIELDSPLIT_SCHUR_PRE_FULL:
190: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from the exact Schur complement\n");
191: break;
192: case PC_FIELDSPLIT_SCHUR_PRE_USER:
193: if (jac->schur_user) {
194: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from user provided matrix\n");
195: } else {
196: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n");
197: }
198: break;
199: default:
200: SETERRQ1(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Invalid Schur preconditioning type: %d", jac->schurpre);
201: }
202: PetscViewerASCIIPrintf(viewer," Split info:\n");
203: PetscViewerASCIIPushTab(viewer);
204: for (i=0; i<jac->nsplits; i++) {
205: if (ilink->fields) {
206: PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);
207: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
208: for (j=0; j<ilink->nfields; j++) {
209: if (j > 0) {
210: PetscViewerASCIIPrintf(viewer,",");
211: }
212: PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);
213: }
214: PetscViewerASCIIPrintf(viewer,"\n");
215: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
216: } else {
217: PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);
218: }
219: ilink = ilink->next;
220: }
221: PetscViewerASCIIPrintf(viewer,"KSP solver for A00 block\n");
222: PetscViewerASCIIPushTab(viewer);
223: if (jac->head) {
224: KSPView(jac->head->ksp,viewer);
225: } else {PetscViewerASCIIPrintf(viewer," not yet available\n");}
226: PetscViewerASCIIPopTab(viewer);
227: if (jac->head && jac->kspupper != jac->head->ksp) {
228: PetscViewerASCIIPrintf(viewer,"KSP solver for upper A00 in upper triangular factor \n");
229: PetscViewerASCIIPushTab(viewer);
230: if (jac->kspupper) {KSPView(jac->kspupper,viewer);}
231: else {PetscViewerASCIIPrintf(viewer," not yet available\n");}
232: PetscViewerASCIIPopTab(viewer);
233: }
234: PetscViewerASCIIPrintf(viewer,"KSP solver for S = A11 - A10 inv(A00) A01 \n");
235: PetscViewerASCIIPushTab(viewer);
236: if (jac->kspschur) {
237: KSPView(jac->kspschur,viewer);
238: } else {
239: PetscViewerASCIIPrintf(viewer," not yet available\n");
240: }
241: PetscViewerASCIIPopTab(viewer);
242: PetscViewerASCIIPopTab(viewer);
243: } else if (isdraw && jac->head) {
244: PetscDraw draw;
245: PetscReal x,y,w,wd,h;
246: PetscInt cnt = 2;
247: char str[32];
249: PetscViewerDrawGetDraw(viewer,0,&draw);
250: PetscDrawGetCurrentPoint(draw,&x,&y);
251: if (jac->kspupper != jac->head->ksp) cnt++;
252: w = 2*PetscMin(1.0 - x,x);
253: wd = w/(cnt + 1);
255: PetscSNPrintf(str,32,"Schur fact. %s",PCFieldSplitSchurFactTypes[jac->schurfactorization]);
256: PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
257: y -= h;
258: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_USER && !jac->schur_user) {
259: PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[PC_FIELDSPLIT_SCHUR_PRE_A11]);
260: } else {
261: PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[jac->schurpre]);
262: }
263: PetscDrawStringBoxed(draw,x+wd*(cnt-1)/2.0,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
264: y -= h;
265: x = x - wd*(cnt-1)/2.0;
267: PetscDrawPushCurrentPoint(draw,x,y);
268: KSPView(jac->head->ksp,viewer);
269: PetscDrawPopCurrentPoint(draw);
270: if (jac->kspupper != jac->head->ksp) {
271: x += wd;
272: PetscDrawPushCurrentPoint(draw,x,y);
273: KSPView(jac->kspupper,viewer);
274: PetscDrawPopCurrentPoint(draw);
275: }
276: x += wd;
277: PetscDrawPushCurrentPoint(draw,x,y);
278: KSPView(jac->kspschur,viewer);
279: PetscDrawPopCurrentPoint(draw);
280: }
281: return(0);
282: }
284: static PetscErrorCode PCView_FieldSplit_GKB(PC pc,PetscViewer viewer)
285: {
286: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
287: PetscErrorCode ierr;
288: PetscBool iascii,isdraw;
289: PetscInt i,j;
290: PC_FieldSplitLink ilink = jac->head;
293: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
294: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
295: if (iascii) {
296: if (jac->bs > 0) {
297: PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs);
298: } else {
299: PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D\n",PCCompositeTypes[jac->type],jac->nsplits);
300: }
301: if (pc->useAmat) {
302: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");
303: }
304: if (jac->diag_use_amat) {
305: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for diagonal blocks\n");
306: }
307: if (jac->offdiag_use_amat) {
308: PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for off-diagonal blocks\n");
309: }
311: PetscViewerASCIIPrintf(viewer," Stopping tolerance=%.1e, delay in error estimate=%D, maximum iterations=%D\n",jac->gkbtol,jac->gkbdelay,jac->gkbmaxit);
312: PetscViewerASCIIPrintf(viewer," Solver info for H = A00 + nu*A01*A01' matrix:\n");
313: PetscViewerASCIIPushTab(viewer);
315: if (ilink->fields) {
316: PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",0);
317: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
318: for (j=0; j<ilink->nfields; j++) {
319: if (j > 0) {
320: PetscViewerASCIIPrintf(viewer,",");
321: }
322: PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);
323: }
324: PetscViewerASCIIPrintf(viewer,"\n");
325: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
326: } else {
327: PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",0);
328: }
329: KSPView(ilink->ksp,viewer);
331: PetscViewerASCIIPopTab(viewer);
332: }
334: if (isdraw) {
335: PetscDraw draw;
336: PetscReal x,y,w,wd;
338: PetscViewerDrawGetDraw(viewer,0,&draw);
339: PetscDrawGetCurrentPoint(draw,&x,&y);
340: w = 2*PetscMin(1.0 - x,x);
341: wd = w/(jac->nsplits + 1);
342: x = x - wd*(jac->nsplits-1)/2.0;
343: for (i=0; i<jac->nsplits; i++) {
344: PetscDrawPushCurrentPoint(draw,x,y);
345: KSPView(ilink->ksp,viewer);
346: PetscDrawPopCurrentPoint(draw);
347: x += wd;
348: ilink = ilink->next;
349: }
350: }
351: return(0);
352: }
355: /* Precondition: jac->bs is set to a meaningful value */
356: static PetscErrorCode PCFieldSplitSetRuntimeSplits_Private(PC pc)
357: {
359: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
360: PetscInt i,nfields,*ifields,nfields_col,*ifields_col;
361: PetscBool flg,flg_col;
362: char optionname[128],splitname[8],optionname_col[128];
365: PetscMalloc1(jac->bs,&ifields);
366: PetscMalloc1(jac->bs,&ifields_col);
367: for (i=0,flg=PETSC_TRUE;; i++) {
368: PetscSNPrintf(splitname,sizeof(splitname),"%D",i);
369: PetscSNPrintf(optionname,sizeof(optionname),"-pc_fieldsplit_%D_fields",i);
370: PetscSNPrintf(optionname_col,sizeof(optionname_col),"-pc_fieldsplit_%D_fields_col",i);
371: nfields = jac->bs;
372: nfields_col = jac->bs;
373: PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,optionname,ifields,&nfields,&flg);
374: PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,optionname_col,ifields_col,&nfields_col,&flg_col);
375: if (!flg) break;
376: else if (flg && !flg_col) {
377: if (!nfields) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields");
378: PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields);
379: } else {
380: if (!nfields || !nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields");
381: if (nfields != nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Number of row and column fields must match");
382: PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields_col);
383: }
384: }
385: if (i > 0) {
386: /* Makes command-line setting of splits take precedence over setting them in code.
387: Otherwise subsequent calls to PCFieldSplitSetIS() or PCFieldSplitSetFields() would
388: create new splits, which would probably not be what the user wanted. */
389: jac->splitdefined = PETSC_TRUE;
390: }
391: PetscFree(ifields);
392: PetscFree(ifields_col);
393: return(0);
394: }
396: static PetscErrorCode PCFieldSplitSetDefaults(PC pc)
397: {
398: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
399: PetscErrorCode ierr;
400: PC_FieldSplitLink ilink = jac->head;
401: PetscBool fieldsplit_default = PETSC_FALSE,coupling = PETSC_FALSE;
402: PetscInt i;
405: /*
406: Kinda messy, but at least this now uses DMCreateFieldDecomposition().
407: Should probably be rewritten.
408: */
409: if (!ilink) {
410: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_coupling",&coupling,NULL);
411: if (pc->dm && jac->dm_splits && !jac->detect && !coupling) {
412: PetscInt numFields, f, i, j;
413: char **fieldNames;
414: IS *fields;
415: DM *dms;
416: DM subdm[128];
417: PetscBool flg;
419: DMCreateFieldDecomposition(pc->dm, &numFields, &fieldNames, &fields, &dms);
420: /* Allow the user to prescribe the splits */
421: for (i = 0, flg = PETSC_TRUE;; i++) {
422: PetscInt ifields[128];
423: IS compField;
424: char optionname[128], splitname[8];
425: PetscInt nfields = numFields;
427: PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%D_fields", i);
428: PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix, optionname, ifields, &nfields, &flg);
429: if (!flg) break;
430: if (numFields > 128) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot currently support %d > 128 fields", numFields);
431: DMCreateSubDM(pc->dm, nfields, ifields, &compField, &subdm[i]);
432: if (nfields == 1) {
433: PCFieldSplitSetIS(pc, fieldNames[ifields[0]], compField);
434: } else {
435: PetscSNPrintf(splitname, sizeof(splitname), "%D", i);
436: PCFieldSplitSetIS(pc, splitname, compField);
437: }
438: ISDestroy(&compField);
439: for (j = 0; j < nfields; ++j) {
440: f = ifields[j];
441: PetscFree(fieldNames[f]);
442: ISDestroy(&fields[f]);
443: }
444: }
445: if (i == 0) {
446: for (f = 0; f < numFields; ++f) {
447: PCFieldSplitSetIS(pc, fieldNames[f], fields[f]);
448: PetscFree(fieldNames[f]);
449: ISDestroy(&fields[f]);
450: }
451: } else {
452: for (j=0; j<numFields; j++) {
453: DMDestroy(dms+j);
454: }
455: PetscFree(dms);
456: PetscMalloc1(i, &dms);
457: for (j = 0; j < i; ++j) dms[j] = subdm[j];
458: }
459: PetscFree(fieldNames);
460: PetscFree(fields);
461: if (dms) {
462: PetscInfo(pc, "Setting up physics based fieldsplit preconditioner using the embedded DM\n");
463: for (ilink = jac->head, i = 0; ilink; ilink = ilink->next, ++i) {
464: const char *prefix;
465: PetscObjectGetOptionsPrefix((PetscObject)(ilink->ksp),&prefix);
466: PetscObjectSetOptionsPrefix((PetscObject)(dms[i]), prefix);
467: KSPSetDM(ilink->ksp, dms[i]);
468: KSPSetDMActive(ilink->ksp, PETSC_FALSE);
469: {
470: PetscErrorCode (*func)(KSP,Mat,Mat,void*);
471: void *ctx;
473: DMKSPGetComputeOperators(pc->dm, &func, &ctx);
474: DMKSPSetComputeOperators(dms[i], func, ctx);
475: }
476: PetscObjectIncrementTabLevel((PetscObject)dms[i],(PetscObject)ilink->ksp,0);
477: DMDestroy(&dms[i]);
478: }
479: PetscFree(dms);
480: }
481: } else {
482: if (jac->bs <= 0) {
483: if (pc->pmat) {
484: MatGetBlockSize(pc->pmat,&jac->bs);
485: } else jac->bs = 1;
486: }
488: if (jac->detect) {
489: IS zerodiags,rest;
490: PetscInt nmin,nmax;
492: MatGetOwnershipRange(pc->mat,&nmin,&nmax);
493: if (jac->diag_use_amat) {
494: MatFindZeroDiagonals(pc->mat,&zerodiags);
495: } else {
496: MatFindZeroDiagonals(pc->pmat,&zerodiags);
497: }
498: ISComplement(zerodiags,nmin,nmax,&rest);
499: PCFieldSplitSetIS(pc,"0",rest);
500: PCFieldSplitSetIS(pc,"1",zerodiags);
501: ISDestroy(&zerodiags);
502: ISDestroy(&rest);
503: } else if (coupling) {
504: IS coupling,rest;
505: PetscInt nmin,nmax;
507: MatGetOwnershipRange(pc->mat,&nmin,&nmax);
508: if (jac->offdiag_use_amat) {
509: MatFindOffBlockDiagonalEntries(pc->mat,&coupling);
510: } else {
511: MatFindOffBlockDiagonalEntries(pc->pmat,&coupling);
512: }
513: ISCreateStride(PetscObjectComm((PetscObject)pc->mat),nmax-nmin,nmin,1,&rest);
514: ISSetIdentity(rest);
515: PCFieldSplitSetIS(pc,"0",rest);
516: PCFieldSplitSetIS(pc,"1",coupling);
517: ISDestroy(&coupling);
518: ISDestroy(&rest);
519: } else {
520: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&fieldsplit_default,NULL);
521: if (!fieldsplit_default) {
522: /* Allow user to set fields from command line, if bs was known at the time of PCSetFromOptions_FieldSplit()
523: then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */
524: PCFieldSplitSetRuntimeSplits_Private(pc);
525: if (jac->splitdefined) {PetscInfo(pc,"Splits defined using the options database\n");}
526: }
527: if ((fieldsplit_default || !jac->splitdefined) && !jac->isrestrict) {
528: Mat M = pc->pmat;
529: PetscBool isnest;
531: PetscInfo(pc,"Using default splitting of fields\n");
532: PetscObjectTypeCompare((PetscObject)pc->pmat,MATNEST,&isnest);
533: if (!isnest) {
534: M = pc->mat;
535: PetscObjectTypeCompare((PetscObject)pc->mat,MATNEST,&isnest);
536: }
537: if (isnest) {
538: IS *fields;
539: PetscInt nf;
541: MatNestGetSize(M,&nf,NULL);
542: PetscMalloc1(nf,&fields);
543: MatNestGetISs(M,fields,NULL);
544: for (i=0;i<nf;i++) {
545: PCFieldSplitSetIS(pc,NULL,fields[i]);
546: }
547: PetscFree(fields);
548: } else {
549: for (i=0; i<jac->bs; i++) {
550: char splitname[8];
551: PetscSNPrintf(splitname,sizeof(splitname),"%D",i);
552: PCFieldSplitSetFields(pc,splitname,1,&i,&i);
553: }
554: jac->defaultsplit = PETSC_TRUE;
555: }
556: }
557: }
558: }
559: } else if (jac->nsplits == 1) {
560: if (ilink->is) {
561: IS is2;
562: PetscInt nmin,nmax;
564: MatGetOwnershipRange(pc->mat,&nmin,&nmax);
565: ISComplement(ilink->is,nmin,nmax,&is2);
566: PCFieldSplitSetIS(pc,"1",is2);
567: ISDestroy(&is2);
568: } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Must provide at least two sets of fields to PCFieldSplit()");
569: }
571: if (jac->nsplits < 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Unhandled case, must have at least two fields, not %d", jac->nsplits);
572: return(0);
573: }
575: static PetscErrorCode MatGolubKahanComputeExplicitOperator(Mat A,Mat B,Mat C,Mat *H,PetscReal gkbnu)
576: {
577: PetscErrorCode ierr;
578: Mat BT,T;
579: PetscReal nrmT,nrmB;
582: MatHermitianTranspose(C,MAT_INITIAL_MATRIX,&T); /* Test if augmented matrix is symmetric */
583: MatAXPY(T,-1.0,B,DIFFERENT_NONZERO_PATTERN);
584: MatNorm(T,NORM_1,&nrmT);
585: MatNorm(B,NORM_1,&nrmB);
586: if (nrmB > 0) {
587: if (nrmT/nrmB >= PETSC_SMALL) {
588: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix is not symmetric/hermitian, GKB is not applicable.");
589: }
590: }
591: /* Compute augmented Lagrangian matrix H = A00 + nu*A01*A01'. This corresponds to */
592: /* setting N := 1/nu*I in [Ar13]. */
593: MatHermitianTranspose(B,MAT_INITIAL_MATRIX,&BT);
594: MatMatMult(B,BT,MAT_INITIAL_MATRIX,PETSC_DEFAULT,H); /* H = A01*A01' */
595: MatAYPX(*H,gkbnu,A,DIFFERENT_NONZERO_PATTERN); /* H = A00 + nu*A01*A01' */
597: MatDestroy(&BT);
598: MatDestroy(&T);
599: return(0);
600: }
602: PETSC_EXTERN PetscErrorCode PetscOptionsFindPairPrefix_Private(PetscOptions,const char pre[], const char name[],const char *value[],PetscBool *flg);
604: static PetscErrorCode PCSetUp_FieldSplit(PC pc)
605: {
606: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
607: PetscErrorCode ierr;
608: PC_FieldSplitLink ilink;
609: PetscInt i,nsplit;
610: PetscBool sorted, sorted_col;
613: pc->failedreason = PC_NOERROR;
614: PCFieldSplitSetDefaults(pc);
615: nsplit = jac->nsplits;
616: ilink = jac->head;
618: /* get the matrices for each split */
619: if (!jac->issetup) {
620: PetscInt rstart,rend,nslots,bs;
622: jac->issetup = PETSC_TRUE;
624: /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */
625: if (jac->defaultsplit || !ilink->is) {
626: if (jac->bs <= 0) jac->bs = nsplit;
627: }
628: bs = jac->bs;
629: MatGetOwnershipRange(pc->pmat,&rstart,&rend);
630: nslots = (rend - rstart)/bs;
631: for (i=0; i<nsplit; i++) {
632: if (jac->defaultsplit) {
633: ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+i,nsplit,&ilink->is);
634: ISDuplicate(ilink->is,&ilink->is_col);
635: } else if (!ilink->is) {
636: if (ilink->nfields > 1) {
637: PetscInt *ii,*jj,j,k,nfields = ilink->nfields,*fields = ilink->fields,*fields_col = ilink->fields_col;
638: PetscMalloc1(ilink->nfields*nslots,&ii);
639: PetscMalloc1(ilink->nfields*nslots,&jj);
640: for (j=0; j<nslots; j++) {
641: for (k=0; k<nfields; k++) {
642: ii[nfields*j + k] = rstart + bs*j + fields[k];
643: jj[nfields*j + k] = rstart + bs*j + fields_col[k];
644: }
645: }
646: ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,ii,PETSC_OWN_POINTER,&ilink->is);
647: ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,jj,PETSC_OWN_POINTER,&ilink->is_col);
648: ISSetBlockSize(ilink->is, nfields);
649: ISSetBlockSize(ilink->is_col, nfields);
650: } else {
651: ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields[0],bs,&ilink->is);
652: ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields_col[0],bs,&ilink->is_col);
653: }
654: }
655: ISSorted(ilink->is,&sorted);
656: if (ilink->is_col) { ISSorted(ilink->is_col,&sorted_col); }
657: if (!sorted || !sorted_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Fields must be sorted when creating split");
658: ilink = ilink->next;
659: }
660: }
662: ilink = jac->head;
663: if (!jac->pmat) {
664: Vec xtmp;
666: MatCreateVecs(pc->pmat,&xtmp,NULL);
667: PetscMalloc1(nsplit,&jac->pmat);
668: PetscMalloc2(nsplit,&jac->x,nsplit,&jac->y);
669: for (i=0; i<nsplit; i++) {
670: MatNullSpace sp;
672: /* Check for preconditioning matrix attached to IS */
673: PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &jac->pmat[i]);
674: if (jac->pmat[i]) {
675: PetscObjectReference((PetscObject) jac->pmat[i]);
676: if (jac->type == PC_COMPOSITE_SCHUR) {
677: jac->schur_user = jac->pmat[i];
679: PetscObjectReference((PetscObject) jac->schur_user);
680: }
681: } else {
682: const char *prefix;
683: MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->pmat[i]);
684: KSPGetOptionsPrefix(ilink->ksp,&prefix);
685: MatSetOptionsPrefix(jac->pmat[i],prefix);
686: MatViewFromOptions(jac->pmat[i],NULL,"-mat_view");
687: }
688: /* create work vectors for each split */
689: MatCreateVecs(jac->pmat[i],&jac->x[i],&jac->y[i]);
690: ilink->x = jac->x[i]; ilink->y = jac->y[i]; ilink->z = NULL;
691: /* compute scatter contexts needed by multiplicative versions and non-default splits */
692: VecScatterCreate(xtmp,ilink->is,jac->x[i],NULL,&ilink->sctx);
693: PetscObjectQuery((PetscObject) ilink->is, "nearnullspace", (PetscObject*) &sp);
694: if (sp) {
695: MatSetNearNullSpace(jac->pmat[i], sp);
696: }
697: ilink = ilink->next;
698: }
699: VecDestroy(&xtmp);
700: } else {
701: MatReuse scall;
702: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
703: for (i=0; i<nsplit; i++) {
704: MatDestroy(&jac->pmat[i]);
705: }
706: scall = MAT_INITIAL_MATRIX;
707: } else scall = MAT_REUSE_MATRIX;
709: for (i=0; i<nsplit; i++) {
710: Mat pmat;
712: /* Check for preconditioning matrix attached to IS */
713: PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &pmat);
714: if (!pmat) {
715: MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is_col,scall,&jac->pmat[i]);
716: }
717: ilink = ilink->next;
718: }
719: }
720: if (jac->diag_use_amat) {
721: ilink = jac->head;
722: if (!jac->mat) {
723: PetscMalloc1(nsplit,&jac->mat);
724: for (i=0; i<nsplit; i++) {
725: MatCreateSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->mat[i]);
726: ilink = ilink->next;
727: }
728: } else {
729: MatReuse scall;
730: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
731: for (i=0; i<nsplit; i++) {
732: MatDestroy(&jac->mat[i]);
733: }
734: scall = MAT_INITIAL_MATRIX;
735: } else scall = MAT_REUSE_MATRIX;
737: for (i=0; i<nsplit; i++) {
738: MatCreateSubMatrix(pc->mat,ilink->is,ilink->is_col,scall,&jac->mat[i]);
739: ilink = ilink->next;
740: }
741: }
742: } else {
743: jac->mat = jac->pmat;
744: }
746: /* Check for null space attached to IS */
747: ilink = jac->head;
748: for (i=0; i<nsplit; i++) {
749: MatNullSpace sp;
751: PetscObjectQuery((PetscObject) ilink->is, "nullspace", (PetscObject*) &sp);
752: if (sp) {
753: MatSetNullSpace(jac->mat[i], sp);
754: }
755: ilink = ilink->next;
756: }
758: if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR && jac->type != PC_COMPOSITE_GKB) {
759: /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */
760: /* FIXME: Can/should we reuse jac->mat whenever (jac->diag_use_amat) is true? */
761: ilink = jac->head;
762: if (nsplit == 2 && jac->type == PC_COMPOSITE_MULTIPLICATIVE) {
763: /* 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 */
764: if (!jac->Afield) {
765: PetscCalloc1(nsplit,&jac->Afield);
766: if (jac->offdiag_use_amat) {
767: MatCreateSubMatrix(pc->mat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);
768: } else {
769: MatCreateSubMatrix(pc->pmat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);
770: }
771: } else {
772: MatReuse scall;
774: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
775: MatDestroy(&jac->Afield[1]);
776: scall = MAT_INITIAL_MATRIX;
777: } else scall = MAT_REUSE_MATRIX;
779: if (jac->offdiag_use_amat) {
780: MatCreateSubMatrix(pc->mat,ilink->next->is,ilink->is,scall,&jac->Afield[1]);
781: } else {
782: MatCreateSubMatrix(pc->pmat,ilink->next->is,ilink->is,scall,&jac->Afield[1]);
783: }
784: }
785: } else {
786: if (!jac->Afield) {
787: PetscMalloc1(nsplit,&jac->Afield);
788: for (i=0; i<nsplit; i++) {
789: if (jac->offdiag_use_amat) {
790: MatCreateSubMatrix(pc->mat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);
791: } else {
792: MatCreateSubMatrix(pc->pmat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);
793: }
794: ilink = ilink->next;
795: }
796: } else {
797: MatReuse scall;
798: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
799: for (i=0; i<nsplit; i++) {
800: MatDestroy(&jac->Afield[i]);
801: }
802: scall = MAT_INITIAL_MATRIX;
803: } else scall = MAT_REUSE_MATRIX;
805: for (i=0; i<nsplit; i++) {
806: if (jac->offdiag_use_amat) {
807: MatCreateSubMatrix(pc->mat,ilink->is,NULL,scall,&jac->Afield[i]);
808: } else {
809: MatCreateSubMatrix(pc->pmat,ilink->is,NULL,scall,&jac->Afield[i]);
810: }
811: ilink = ilink->next;
812: }
813: }
814: }
815: }
817: if (jac->type == PC_COMPOSITE_SCHUR) {
818: IS ccis;
819: PetscBool isspd;
820: PetscInt rstart,rend;
821: char lscname[256];
822: PetscObject LSC_L;
824: if (nsplit != 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"To use Schur complement preconditioner you must have exactly 2 fields");
826: /* If pc->mat is SPD, don't scale by -1 the Schur complement */
827: if (jac->schurscale == (PetscScalar)-1.0) {
828: MatGetOption(pc->pmat,MAT_SPD,&isspd);
829: jac->schurscale = (isspd == PETSC_TRUE) ? 1.0 : -1.0;
830: }
832: /* When extracting off-diagonal submatrices, we take complements from this range */
833: MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);
835: if (jac->schur) {
836: KSP kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper;
837: MatReuse scall;
839: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
840: scall = MAT_INITIAL_MATRIX;
841: MatDestroy(&jac->B);
842: MatDestroy(&jac->C);
843: } else scall = MAT_REUSE_MATRIX;
845: MatSchurComplementGetKSP(jac->schur, &kspInner);
846: ilink = jac->head;
847: ISComplement(ilink->is_col,rstart,rend,&ccis);
848: if (jac->offdiag_use_amat) {
849: MatCreateSubMatrix(pc->mat,ilink->is,ccis,scall,&jac->B);
850: } else {
851: MatCreateSubMatrix(pc->pmat,ilink->is,ccis,scall,&jac->B);
852: }
853: ISDestroy(&ccis);
854: ilink = ilink->next;
855: ISComplement(ilink->is_col,rstart,rend,&ccis);
856: if (jac->offdiag_use_amat) {
857: MatCreateSubMatrix(pc->mat,ilink->is,ccis,scall,&jac->C);
858: } else {
859: MatCreateSubMatrix(pc->pmat,ilink->is,ccis,scall,&jac->C);
860: }
861: ISDestroy(&ccis);
862: MatSchurComplementUpdateSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);
863: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
864: MatDestroy(&jac->schurp);
865: MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);
866: }
867: if (kspA != kspInner) {
868: KSPSetOperators(kspA,jac->mat[0],jac->pmat[0]);
869: }
870: if (kspUpper != kspA) {
871: KSPSetOperators(kspUpper,jac->mat[0],jac->pmat[0]);
872: }
873: KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));
874: } else {
875: const char *Dprefix;
876: char schurprefix[256], schurmatprefix[256];
877: char schurtestoption[256];
878: MatNullSpace sp;
879: PetscBool flg;
880: KSP kspt;
882: /* extract the A01 and A10 matrices */
883: ilink = jac->head;
884: ISComplement(ilink->is_col,rstart,rend,&ccis);
885: if (jac->offdiag_use_amat) {
886: MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
887: } else {
888: MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
889: }
890: ISDestroy(&ccis);
891: ilink = ilink->next;
892: ISComplement(ilink->is_col,rstart,rend,&ccis);
893: if (jac->offdiag_use_amat) {
894: MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
895: } else {
896: MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
897: }
898: ISDestroy(&ccis);
900: /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */
901: MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur);
902: MatSetType(jac->schur,MATSCHURCOMPLEMENT);
903: MatSchurComplementSetSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);
904: PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
905: MatSetOptionsPrefix(jac->schur,schurmatprefix);
906: MatSchurComplementGetKSP(jac->schur,&kspt);
907: KSPSetOptionsPrefix(kspt,schurmatprefix);
909: /* Note: this is not true in general */
910: MatGetNullSpace(jac->mat[1], &sp);
911: if (sp) {
912: MatSetNullSpace(jac->schur, sp);
913: }
915: PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);
916: PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
917: if (flg) {
918: DM dmInner;
919: KSP kspInner;
920: PC pcInner;
922: MatSchurComplementGetKSP(jac->schur, &kspInner);
923: KSPReset(kspInner);
924: KSPSetOperators(kspInner,jac->mat[0],jac->pmat[0]);
925: PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
926: /* Indent this deeper to emphasize the "inner" nature of this solver. */
927: PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject) pc, 2);
928: PetscObjectIncrementTabLevel((PetscObject)kspInner->pc, (PetscObject) pc, 2);
929: KSPSetOptionsPrefix(kspInner, schurprefix);
931: /* Set DM for new solver */
932: KSPGetDM(jac->head->ksp, &dmInner);
933: KSPSetDM(kspInner, dmInner);
934: KSPSetDMActive(kspInner, PETSC_FALSE);
936: /* Defaults to PCKSP as preconditioner */
937: KSPGetPC(kspInner, &pcInner);
938: PCSetType(pcInner, PCKSP);
939: PCKSPSetKSP(pcInner, jac->head->ksp);
940: } else {
941: /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or
942: * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact,
943: * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for
944: * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make
945: * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used
946: * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */
947: KSPSetType(jac->head->ksp,KSPGMRES);
948: MatSchurComplementSetKSP(jac->schur,jac->head->ksp);
949: }
950: KSPSetOperators(jac->head->ksp,jac->mat[0],jac->pmat[0]);
951: KSPSetFromOptions(jac->head->ksp);
952: MatSetFromOptions(jac->schur);
954: PetscObjectTypeCompare((PetscObject)jac->schur, MATSCHURCOMPLEMENT, &flg);
955: if (flg) { /* Need to do this otherwise PCSetUp_KSP will overwrite the amat of jac->head->ksp */
956: KSP kspInner;
957: PC pcInner;
959: MatSchurComplementGetKSP(jac->schur, &kspInner);
960: KSPGetPC(kspInner, &pcInner);
961: PetscObjectTypeCompare((PetscObject)pcInner, PCKSP, &flg);
962: if (flg) {
963: KSP ksp;
965: PCKSPGetKSP(pcInner, &ksp);
966: if (ksp == jac->head->ksp) {
967: PCSetUseAmat(pcInner, PETSC_TRUE);
968: }
969: }
970: }
971: PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);
972: PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
973: if (flg) {
974: DM dmInner;
976: PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
977: KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper);
978: KSPSetErrorIfNotConverged(jac->kspupper,pc->erroriffailure);
979: KSPSetOptionsPrefix(jac->kspupper, schurprefix);
980: PetscObjectIncrementTabLevel((PetscObject)jac->kspupper, (PetscObject) pc, 1);
981: PetscObjectIncrementTabLevel((PetscObject)jac->kspupper->pc, (PetscObject) pc, 1);
982: KSPGetDM(jac->head->ksp, &dmInner);
983: KSPSetDM(jac->kspupper, dmInner);
984: KSPSetDMActive(jac->kspupper, PETSC_FALSE);
985: KSPSetFromOptions(jac->kspupper);
986: KSPSetOperators(jac->kspupper,jac->mat[0],jac->pmat[0]);
987: VecDuplicate(jac->head->x, &jac->head->z);
988: } else {
989: jac->kspupper = jac->head->ksp;
990: PetscObjectReference((PetscObject) jac->head->ksp);
991: }
993: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
994: MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);
995: }
996: KSPCreate(PetscObjectComm((PetscObject)pc),&jac->kspschur);
997: KSPSetErrorIfNotConverged(jac->kspschur,pc->erroriffailure);
998: PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);
999: PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);
1000: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) {
1001: PC pcschur;
1002: KSPGetPC(jac->kspschur,&pcschur);
1003: PCSetType(pcschur,PCNONE);
1004: /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */
1005: } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) {
1006: MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user);
1007: }
1008: KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));
1009: KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix);
1010: KSPSetOptionsPrefix(jac->kspschur, Dprefix);
1011: /* propagate DM */
1012: {
1013: DM sdm;
1014: KSPGetDM(jac->head->next->ksp, &sdm);
1015: if (sdm) {
1016: KSPSetDM(jac->kspschur, sdm);
1017: KSPSetDMActive(jac->kspschur, PETSC_FALSE);
1018: }
1019: }
1020: /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
1021: /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */
1022: KSPSetFromOptions(jac->kspschur);
1023: }
1024: MatAssemblyBegin(jac->schur,MAT_FINAL_ASSEMBLY);
1025: MatAssemblyEnd(jac->schur,MAT_FINAL_ASSEMBLY);
1027: /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */
1028: PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_L",ilink->splitname);
1029: PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);
1030: if (!LSC_L) {PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);}
1031: if (LSC_L) {PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L);}
1032: PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_Lp",ilink->splitname);
1033: PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);
1034: if (!LSC_L) {PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);}
1035: if (LSC_L) {PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L);}
1036: } else if (jac->type == PC_COMPOSITE_GKB) {
1037: IS ccis;
1038: PetscInt rstart,rend;
1040: if (nsplit != 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"To use GKB preconditioner you must have exactly 2 fields");
1042: ilink = jac->head;
1044: /* When extracting off-diagonal submatrices, we take complements from this range */
1045: MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);
1047: ISComplement(ilink->is_col,rstart,rend,&ccis);
1048: if (jac->offdiag_use_amat) {
1049: MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
1050: } else {
1051: MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
1052: }
1053: ISDestroy(&ccis);
1054: /* Create work vectors for GKB algorithm */
1055: VecDuplicate(ilink->x,&jac->u);
1056: VecDuplicate(ilink->x,&jac->Hu);
1057: VecDuplicate(ilink->x,&jac->w2);
1058: ilink = ilink->next;
1059: ISComplement(ilink->is_col,rstart,rend,&ccis);
1060: if (jac->offdiag_use_amat) {
1061: MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
1062: } else {
1063: MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
1064: }
1065: ISDestroy(&ccis);
1066: /* Create work vectors for GKB algorithm */
1067: VecDuplicate(ilink->x,&jac->v);
1068: VecDuplicate(ilink->x,&jac->d);
1069: VecDuplicate(ilink->x,&jac->w1);
1070: MatGolubKahanComputeExplicitOperator(jac->mat[0],jac->B,jac->C,&jac->H,jac->gkbnu);
1071: PetscCalloc1(jac->gkbdelay,&jac->vecz);
1073: ilink = jac->head;
1074: KSPSetOperators(ilink->ksp,jac->H,jac->H);
1075: if (!jac->suboptionsset) {KSPSetFromOptions(ilink->ksp);}
1076: /* Create gkb_monitor context */
1077: if (jac->gkbmonitor) {
1078: PetscInt tablevel;
1079: PetscViewerCreate(PETSC_COMM_WORLD,&jac->gkbviewer);
1080: PetscViewerSetType(jac->gkbviewer,PETSCVIEWERASCII);
1081: PetscObjectGetTabLevel((PetscObject)ilink->ksp,&tablevel);
1082: PetscViewerASCIISetTab(jac->gkbviewer,tablevel);
1083: PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)ilink->ksp,1);
1084: }
1085: } else {
1086: /* set up the individual splits' PCs */
1087: i = 0;
1088: ilink = jac->head;
1089: while (ilink) {
1090: KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i]);
1091: /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
1092: if (!jac->suboptionsset) {KSPSetFromOptions(ilink->ksp);}
1093: i++;
1094: ilink = ilink->next;
1095: }
1096: }
1098: jac->suboptionsset = PETSC_TRUE;
1099: return(0);
1100: }
1102: #define FieldSplitSplitSolveAdd(ilink,xx,yy) \
1103: (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
1104: VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
1105: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\
1106: KSPSolve(ilink->ksp,ilink->x,ilink->y) || \
1107: KSPCheckSolve(ilink->ksp,pc,ilink->y) || \
1108: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\
1109: VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \
1110: VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE))
1112: static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y)
1113: {
1114: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1115: PetscErrorCode ierr;
1116: PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
1117: KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;
1120: switch (jac->schurfactorization) {
1121: case PC_FIELDSPLIT_SCHUR_FACT_DIAG:
1122: /* [A00 0; 0 -S], positive definite, suitable for MINRES */
1123: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1124: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
1125: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1126: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1127: KSPSolve(kspA,ilinkA->x,ilinkA->y);
1128: KSPCheckSolve(kspA,pc,ilinkA->y);
1129: PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1130: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1131: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
1132: PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1133: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
1134: KSPCheckSolve(jac->kspschur,pc,ilinkD->y);
1135: PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1136: VecScale(ilinkD->y,jac->schurscale);
1137: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1138: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1139: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1140: break;
1141: case PC_FIELDSPLIT_SCHUR_FACT_LOWER:
1142: /* [A00 0; A10 S], suitable for left preconditioning */
1143: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1144: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1145: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1146: KSPSolve(kspA,ilinkA->x,ilinkA->y);
1147: KSPCheckSolve(kspA,pc,ilinkA->y);
1148: PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1149: MatMult(jac->C,ilinkA->y,ilinkD->x);
1150: VecScale(ilinkD->x,-1.);
1151: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
1152: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1153: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
1154: PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1155: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
1156: KSPCheckSolve(jac->kspschur,pc,ilinkD->y);
1157: PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1158: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1159: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1160: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1161: break;
1162: case PC_FIELDSPLIT_SCHUR_FACT_UPPER:
1163: /* [A00 A01; 0 S], suitable for right preconditioning */
1164: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
1165: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
1166: PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1167: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
1168: KSPCheckSolve(jac->kspschur,pc,ilinkD->y);
1169: PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL); MatMult(jac->B,ilinkD->y,ilinkA->x);
1170: VecScale(ilinkA->x,-1.);
1171: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);
1172: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1173: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);
1174: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1175: KSPSolve(kspA,ilinkA->x,ilinkA->y);
1176: KSPCheckSolve(kspA,pc,ilinkA->y);
1177: PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1178: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1179: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1180: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1181: break;
1182: case PC_FIELDSPLIT_SCHUR_FACT_FULL:
1183: /* [1 0; A10 A00^{-1} 1] [A00 0; 0 S] [1 A00^{-1}A01; 0 1] */
1184: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1185: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1186: PetscLogEventBegin(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);
1187: KSPSolve(kspLower,ilinkA->x,ilinkA->y);
1188: KSPCheckSolve(kspLower,pc,ilinkA->y);
1189: PetscLogEventEnd(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);
1190: MatMult(jac->C,ilinkA->y,ilinkD->x);
1191: VecScale(ilinkD->x,-1.0);
1192: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
1193: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
1195: PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1196: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
1197: KSPCheckSolve(jac->kspschur,pc,ilinkD->y);
1198: PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1199: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1201: if (kspUpper == kspA) {
1202: MatMult(jac->B,ilinkD->y,ilinkA->y);
1203: VecAXPY(ilinkA->x,-1.0,ilinkA->y);
1204: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1205: KSPSolve(kspA,ilinkA->x,ilinkA->y);
1206: KSPCheckSolve(kspA,pc,ilinkA->y);
1207: PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1208: } else {
1209: PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1210: KSPSolve(kspA,ilinkA->x,ilinkA->y);
1211: KSPCheckSolve(kspA,pc,ilinkA->y);
1212: MatMult(jac->B,ilinkD->y,ilinkA->x);
1213: PetscLogEventBegin(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);
1214: KSPSolve(kspUpper,ilinkA->x,ilinkA->z);
1215: KSPCheckSolve(kspUpper,pc,ilinkA->z);
1216: PetscLogEventEnd(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);
1217: VecAXPY(ilinkA->y,-1.0,ilinkA->z);
1218: }
1219: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1220: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1221: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1222: }
1223: return(0);
1224: }
1226: static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y)
1227: {
1228: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1229: PetscErrorCode ierr;
1230: PC_FieldSplitLink ilink = jac->head;
1231: PetscInt cnt,bs;
1234: if (jac->type == PC_COMPOSITE_ADDITIVE) {
1235: if (jac->defaultsplit) {
1236: VecGetBlockSize(x,&bs);
1237: 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);
1238: VecGetBlockSize(y,&bs);
1239: 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);
1240: VecStrideGatherAll(x,jac->x,INSERT_VALUES);
1241: while (ilink) {
1242: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1243: KSPSolve(ilink->ksp,ilink->x,ilink->y);
1244: KSPCheckSolve(ilink->ksp,pc,ilink->y);
1245: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1246: ilink = ilink->next;
1247: }
1248: VecStrideScatterAll(jac->y,y,INSERT_VALUES);
1249: } else {
1250: VecSet(y,0.0);
1251: while (ilink) {
1252: FieldSplitSplitSolveAdd(ilink,x,y);
1253: ilink = ilink->next;
1254: }
1255: }
1256: } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) {
1257: VecSet(y,0.0);
1258: /* solve on first block for first block variables */
1259: VecScatterBegin(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);
1260: VecScatterEnd(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);
1261: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1262: KSPSolve(ilink->ksp,ilink->x,ilink->y);
1263: KSPCheckSolve(ilink->ksp,pc,ilink->y);
1264: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1265: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1266: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1268: /* compute the residual only onto second block variables using first block variables */
1269: MatMult(jac->Afield[1],ilink->y,ilink->next->x);
1270: ilink = ilink->next;
1271: VecScale(ilink->x,-1.0);
1272: VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1273: VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1275: /* solve on second block variables */
1276: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1277: KSPSolve(ilink->ksp,ilink->x,ilink->y);
1278: KSPCheckSolve(ilink->ksp,pc,ilink->y);
1279: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1280: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1281: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1282: } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1283: if (!jac->w1) {
1284: VecDuplicate(x,&jac->w1);
1285: VecDuplicate(x,&jac->w2);
1286: }
1287: VecSet(y,0.0);
1288: FieldSplitSplitSolveAdd(ilink,x,y);
1289: cnt = 1;
1290: while (ilink->next) {
1291: ilink = ilink->next;
1292: /* compute the residual only over the part of the vector needed */
1293: MatMult(jac->Afield[cnt++],y,ilink->x);
1294: VecScale(ilink->x,-1.0);
1295: VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1296: VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1297: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1298: KSPSolve(ilink->ksp,ilink->x,ilink->y);
1299: KSPCheckSolve(ilink->ksp,pc,ilink->y);
1300: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1301: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1302: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1303: }
1304: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1305: cnt -= 2;
1306: while (ilink->previous) {
1307: ilink = ilink->previous;
1308: /* compute the residual only over the part of the vector needed */
1309: MatMult(jac->Afield[cnt--],y,ilink->x);
1310: VecScale(ilink->x,-1.0);
1311: VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1312: VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1313: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1314: KSPSolve(ilink->ksp,ilink->x,ilink->y);
1315: KSPCheckSolve(ilink->ksp,pc,ilink->y);
1316: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1317: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1318: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1319: }
1320: }
1321: } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type);
1322: return(0);
1323: }
1326: static PetscErrorCode PCApply_FieldSplit_GKB(PC pc,Vec x,Vec y)
1327: {
1328: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1329: PetscErrorCode ierr;
1330: PC_FieldSplitLink ilinkA = jac->head,ilinkD = ilinkA->next;
1331: KSP ksp = ilinkA->ksp;
1332: Vec u,v,Hu,d,work1,work2;
1333: PetscScalar alpha,z,nrmz2,*vecz;
1334: PetscReal lowbnd,nu,beta;
1335: PetscInt j,iterGKB;
1338: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1339: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
1340: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1341: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
1343: u = jac->u;
1344: v = jac->v;
1345: Hu = jac->Hu;
1346: d = jac->d;
1347: work1 = jac->w1;
1348: work2 = jac->w2;
1349: vecz = jac->vecz;
1351: /* Change RHS to comply with matrix regularization H = A + nu*B*B' */
1352: /* Add q = q + nu*B*b */
1353: if (jac->gkbnu) {
1354: nu = jac->gkbnu;
1355: VecScale(ilinkD->x,jac->gkbnu);
1356: MatMultAdd(jac->B,ilinkD->x,ilinkA->x,ilinkA->x); /* q = q + nu*B*b */
1357: } else {
1358: /* Situation when no augmented Lagrangian is used. Then we set inner */
1359: /* matrix N = I in [Ar13], and thus nu = 1. */
1360: nu = 1;
1361: }
1363: /* Transform rhs from [q,tilde{b}] to [0,b] */
1364: PetscLogEventBegin(ilinkA->event,ksp,ilinkA->x,ilinkA->y,NULL);
1365: KSPSolve(ksp,ilinkA->x,ilinkA->y);
1366: KSPCheckSolve(ksp,pc,ilinkA->y);
1367: PetscLogEventEnd(ilinkA->event,ksp,ilinkA->x,ilinkA->y,NULL);
1368: MatMultHermitianTranspose(jac->B,ilinkA->y,work1);
1369: VecAXPBY(work1,1.0/nu,-1.0,ilinkD->x); /* c = b - B'*x */
1371: /* First step of algorithm */
1372: VecNorm(work1,NORM_2,&beta); /* beta = sqrt(nu*c'*c)*/
1373: KSPCheckDot(ksp,beta);
1374: beta = PetscSqrtScalar(nu)*beta;
1375: VecAXPBY(v,nu/beta,0.0,work1); /* v = nu/beta *c */
1376: MatMult(jac->B,v,work2); /* u = H^{-1}*B*v */
1377: PetscLogEventBegin(ilinkA->event,ksp,work2,u,NULL);
1378: KSPSolve(ksp,work2,u);
1379: KSPCheckSolve(ksp,pc,u);
1380: PetscLogEventEnd(ilinkA->event,ksp,work2,u,NULL);
1381: MatMult(jac->H,u,Hu); /* alpha = u'*H*u */
1382: VecDot(Hu,u,&alpha);
1383: KSPCheckDot(ksp,alpha);
1384: if (PetscRealPart(alpha) <= 0.0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_NOT_CONVERGED,"GKB preconditioner diverged, H is not positive definite");
1385: alpha = PetscSqrtScalar(PetscAbsScalar(alpha));
1386: VecScale(u,1.0/alpha);
1387: VecAXPBY(d,1.0/alpha,0.0,v); /* v = nu/beta *c */
1389: z = beta/alpha;
1390: vecz[1] = z;
1392: /* Computation of first iterate x(1) and p(1) */
1393: VecAXPY(ilinkA->y,z,u);
1394: VecCopy(d,ilinkD->y);
1395: VecScale(ilinkD->y,-z);
1397: iterGKB = 1; lowbnd = 2*jac->gkbtol;
1398: if (jac->gkbmonitor) {
1399: PetscViewerASCIIPrintf(jac->gkbviewer,"%3D GKB Lower bound estimate %14.12e\n",iterGKB,lowbnd);
1400: }
1402: while (iterGKB < jac->gkbmaxit && lowbnd > jac->gkbtol) {
1403: iterGKB += 1;
1404: MatMultHermitianTranspose(jac->B,u,work1); /* v <- nu*(B'*u-alpha/nu*v) */
1405: VecAXPBY(v,nu,-alpha,work1);
1406: VecNorm(v,NORM_2,&beta); /* beta = sqrt(nu)*v'*v */
1407: beta = beta/PetscSqrtScalar(nu);
1408: VecScale(v,1.0/beta);
1409: MatMult(jac->B,v,work2); /* u <- H^{-1}*(B*v-beta*H*u) */
1410: MatMult(jac->H,u,Hu);
1411: VecAXPY(work2,-beta,Hu);
1412: PetscLogEventBegin(ilinkA->event,ksp,work2,u,NULL);
1413: KSPSolve(ksp,work2,u);
1414: KSPCheckSolve(ksp,pc,u);
1415: PetscLogEventEnd(ilinkA->event,ksp,work2,u,NULL);
1416: MatMult(jac->H,u,Hu); /* alpha = u'*H*u */
1417: VecDot(Hu,u,&alpha);
1418: KSPCheckDot(ksp,alpha);
1419: if (PetscRealPart(alpha) <= 0.0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_NOT_CONVERGED,"GKB preconditioner diverged, H is not positive definite");
1420: alpha = PetscSqrtScalar(PetscAbsScalar(alpha));
1421: VecScale(u,1.0/alpha);
1423: z = -beta/alpha*z; /* z <- beta/alpha*z */
1424: vecz[0] = z;
1426: /* Computation of new iterate x(i+1) and p(i+1) */
1427: VecAXPBY(d,1.0/alpha,-beta/alpha,v); /* d = (v-beta*d)/alpha */
1428: VecAXPY(ilinkA->y,z,u); /* r = r + z*u */
1429: VecAXPY(ilinkD->y,-z,d); /* p = p - z*d */
1430: MatMult(jac->H,ilinkA->y,Hu); /* ||u||_H = u'*H*u */
1431: VecDot(Hu,ilinkA->y,&nrmz2);
1433: /* Compute Lower Bound estimate */
1434: if (iterGKB > jac->gkbdelay) {
1435: lowbnd = 0.0;
1436: for (j=0; j<jac->gkbdelay; j++) {
1437: lowbnd += PetscAbsScalar(vecz[j]*vecz[j]);
1438: }
1439: lowbnd = PetscSqrtScalar(lowbnd/PetscAbsScalar(nrmz2));
1440: }
1442: for (j=0; j<jac->gkbdelay-1; j++) {
1443: vecz[jac->gkbdelay-j-1] = vecz[jac->gkbdelay-j-2];
1444: }
1445: if (jac->gkbmonitor) {
1446: PetscViewerASCIIPrintf(jac->gkbviewer,"%3D GKB Lower bound estimate %14.12e\n",iterGKB,lowbnd);
1447: }
1448: }
1450: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1451: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1452: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1453: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1455: return(0);
1456: }
1459: #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \
1460: (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
1461: VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
1462: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \
1463: KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \
1464: KSPCheckSolve(ilink->ksp,pc,ilink->x) || \
1465: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \
1466: VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \
1467: VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE))
1469: static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y)
1470: {
1471: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1472: PetscErrorCode ierr;
1473: PC_FieldSplitLink ilink = jac->head;
1474: PetscInt bs;
1477: if (jac->type == PC_COMPOSITE_ADDITIVE) {
1478: if (jac->defaultsplit) {
1479: VecGetBlockSize(x,&bs);
1480: 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);
1481: VecGetBlockSize(y,&bs);
1482: 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);
1483: VecStrideGatherAll(x,jac->x,INSERT_VALUES);
1484: while (ilink) {
1485: PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1486: KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);
1487: KSPCheckSolve(ilink->ksp,pc,ilink->y);
1488: PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1489: ilink = ilink->next;
1490: }
1491: VecStrideScatterAll(jac->y,y,INSERT_VALUES);
1492: } else {
1493: VecSet(y,0.0);
1494: while (ilink) {
1495: FieldSplitSplitSolveAddTranspose(ilink,x,y);
1496: ilink = ilink->next;
1497: }
1498: }
1499: } else {
1500: if (!jac->w1) {
1501: VecDuplicate(x,&jac->w1);
1502: VecDuplicate(x,&jac->w2);
1503: }
1504: VecSet(y,0.0);
1505: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1506: FieldSplitSplitSolveAddTranspose(ilink,x,y);
1507: while (ilink->next) {
1508: ilink = ilink->next;
1509: MatMultTranspose(pc->mat,y,jac->w1);
1510: VecWAXPY(jac->w2,-1.0,jac->w1,x);
1511: FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1512: }
1513: while (ilink->previous) {
1514: ilink = ilink->previous;
1515: MatMultTranspose(pc->mat,y,jac->w1);
1516: VecWAXPY(jac->w2,-1.0,jac->w1,x);
1517: FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1518: }
1519: } else {
1520: while (ilink->next) { /* get to last entry in linked list */
1521: ilink = ilink->next;
1522: }
1523: FieldSplitSplitSolveAddTranspose(ilink,x,y);
1524: while (ilink->previous) {
1525: ilink = ilink->previous;
1526: MatMultTranspose(pc->mat,y,jac->w1);
1527: VecWAXPY(jac->w2,-1.0,jac->w1,x);
1528: FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1529: }
1530: }
1531: }
1532: return(0);
1533: }
1535: static PetscErrorCode PCReset_FieldSplit(PC pc)
1536: {
1537: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1538: PetscErrorCode ierr;
1539: PC_FieldSplitLink ilink = jac->head,next;
1542: while (ilink) {
1543: KSPDestroy(&ilink->ksp);
1544: VecDestroy(&ilink->x);
1545: VecDestroy(&ilink->y);
1546: VecDestroy(&ilink->z);
1547: VecScatterDestroy(&ilink->sctx);
1548: ISDestroy(&ilink->is);
1549: ISDestroy(&ilink->is_col);
1550: PetscFree(ilink->splitname);
1551: PetscFree(ilink->fields);
1552: PetscFree(ilink->fields_col);
1553: next = ilink->next;
1554: PetscFree(ilink);
1555: ilink = next;
1556: }
1557: jac->head = NULL;
1558: PetscFree2(jac->x,jac->y);
1559: if (jac->mat && jac->mat != jac->pmat) {
1560: MatDestroyMatrices(jac->nsplits,&jac->mat);
1561: } else if (jac->mat) {
1562: jac->mat = NULL;
1563: }
1564: if (jac->pmat) {MatDestroyMatrices(jac->nsplits,&jac->pmat);}
1565: if (jac->Afield) {MatDestroyMatrices(jac->nsplits,&jac->Afield);}
1566: jac->nsplits = 0;
1567: VecDestroy(&jac->w1);
1568: VecDestroy(&jac->w2);
1569: MatDestroy(&jac->schur);
1570: MatDestroy(&jac->schurp);
1571: MatDestroy(&jac->schur_user);
1572: KSPDestroy(&jac->kspschur);
1573: KSPDestroy(&jac->kspupper);
1574: MatDestroy(&jac->B);
1575: MatDestroy(&jac->C);
1576: MatDestroy(&jac->H);
1577: VecDestroy(&jac->u);
1578: VecDestroy(&jac->v);
1579: VecDestroy(&jac->Hu);
1580: VecDestroy(&jac->d);
1581: PetscFree(jac->vecz);
1582: PetscViewerDestroy(&jac->gkbviewer);
1583: jac->isrestrict = PETSC_FALSE;
1584: return(0);
1585: }
1587: static PetscErrorCode PCDestroy_FieldSplit(PC pc)
1588: {
1589: PetscErrorCode ierr;
1592: PCReset_FieldSplit(pc);
1593: PetscFree(pc->data);
1594: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurGetSubKSP_C",NULL);
1595: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL);
1596: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL);
1597: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL);
1598: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL);
1599: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL);
1600: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL);
1601: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL);
1602: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL);
1603: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",NULL);
1604: return(0);
1605: }
1607: static PetscErrorCode PCSetFromOptions_FieldSplit(PetscOptionItems *PetscOptionsObject,PC pc)
1608: {
1609: PetscErrorCode ierr;
1610: PetscInt bs;
1611: PetscBool flg;
1612: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1613: PCCompositeType ctype;
1616: PetscOptionsHead(PetscOptionsObject,"FieldSplit options");
1617: PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL);
1618: PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);
1619: if (flg) {
1620: PCFieldSplitSetBlockSize(pc,bs);
1621: }
1622: jac->diag_use_amat = pc->useAmat;
1623: 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);
1624: jac->offdiag_use_amat = pc->useAmat;
1625: 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);
1626: PetscOptionsBool("-pc_fieldsplit_detect_saddle_point","Form 2-way split by detecting zero diagonal entries", "PCFieldSplitSetDetectSaddlePoint",jac->detect,&jac->detect,NULL);
1627: PCFieldSplitSetDetectSaddlePoint(pc,jac->detect); /* Sets split type and Schur PC type */
1628: PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);
1629: if (flg) {
1630: PCFieldSplitSetType(pc,ctype);
1631: }
1632: /* Only setup fields once */
1633: if ((jac->bs > 0) && (jac->nsplits == 0)) {
1634: /* only allow user to set fields from command line if bs is already known.
1635: otherwise user can set them in PCFieldSplitSetDefaults() */
1636: PCFieldSplitSetRuntimeSplits_Private(pc);
1637: if (jac->splitdefined) {PetscInfo(pc,"Splits defined using the options database\n");}
1638: }
1639: if (jac->type == PC_COMPOSITE_SCHUR) {
1640: PetscOptionsGetEnum(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);
1641: if (flg) {PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");}
1642: 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);
1643: PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSetSchurPre",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL);
1644: PetscOptionsScalar("-pc_fieldsplit_schur_scale","Scale Schur complement","PCFieldSplitSetSchurScale",jac->schurscale,&jac->schurscale,NULL);
1645: } else if (jac->type == PC_COMPOSITE_GKB) {
1646: PetscOptionsReal("-pc_fieldsplit_gkb_tol","The tolerance for the lower bound stopping criterion","PCFieldSplitGKBTol",jac->gkbtol,&jac->gkbtol,NULL);
1647: PetscOptionsInt("-pc_fieldsplit_gkb_delay","The delay value for lower bound criterion","PCFieldSplitGKBDelay",jac->gkbdelay,&jac->gkbdelay,NULL);
1648: PetscOptionsReal("-pc_fieldsplit_gkb_nu","Parameter in augmented Lagrangian approach","PCFieldSplitGKBNu",jac->gkbnu,&jac->gkbnu,NULL);
1649: if (jac->gkbnu < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nu cannot be less than 0: value %f",jac->gkbnu);
1650: PetscOptionsInt("-pc_fieldsplit_gkb_maxit","Maximum allowed number of iterations","PCFieldSplitGKBMaxit",jac->gkbmaxit,&jac->gkbmaxit,NULL);
1651: PetscOptionsBool("-pc_fieldsplit_gkb_monitor","Prints number of GKB iterations and error","PCFieldSplitGKB",jac->gkbmonitor,&jac->gkbmonitor,NULL);
1652: }
1653: PetscOptionsTail();
1654: return(0);
1655: }
1657: /*------------------------------------------------------------------------------------*/
1659: static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1660: {
1661: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1662: PetscErrorCode ierr;
1663: PC_FieldSplitLink ilink,next = jac->head;
1664: char prefix[128];
1665: PetscInt i;
1668: if (jac->splitdefined) {
1669: PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1670: return(0);
1671: }
1672: for (i=0; i<n; i++) {
1673: if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs);
1674: if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]);
1675: }
1676: PetscNew(&ilink);
1677: if (splitname) {
1678: PetscStrallocpy(splitname,&ilink->splitname);
1679: } else {
1680: PetscMalloc1(3,&ilink->splitname);
1681: PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);
1682: }
1683: 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 */
1684: PetscMalloc1(n,&ilink->fields);
1685: PetscArraycpy(ilink->fields,fields,n);
1686: PetscMalloc1(n,&ilink->fields_col);
1687: PetscArraycpy(ilink->fields_col,fields_col,n);
1689: ilink->nfields = n;
1690: ilink->next = NULL;
1691: KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1692: KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);
1693: PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1694: KSPSetType(ilink->ksp,KSPPREONLY);
1695: PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);
1697: PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);
1698: KSPSetOptionsPrefix(ilink->ksp,prefix);
1700: if (!next) {
1701: jac->head = ilink;
1702: ilink->previous = NULL;
1703: } else {
1704: while (next->next) {
1705: next = next->next;
1706: }
1707: next->next = ilink;
1708: ilink->previous = next;
1709: }
1710: jac->nsplits++;
1711: return(0);
1712: }
1714: static PetscErrorCode PCFieldSplitSchurGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp)
1715: {
1716: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1720: *subksp = NULL;
1721: if (n) *n = 0;
1722: if (jac->type == PC_COMPOSITE_SCHUR) {
1723: PetscInt nn;
1725: if (!jac->schur) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitSchurGetSubKSP()");
1726: if (jac->nsplits != 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Unexpected number of splits %D != 2",jac->nsplits);
1727: nn = jac->nsplits + (jac->kspupper != jac->head->ksp ? 1 : 0);
1728: PetscMalloc1(nn,subksp);
1729: (*subksp)[0] = jac->head->ksp;
1730: (*subksp)[1] = jac->kspschur;
1731: if (jac->kspupper != jac->head->ksp) (*subksp)[2] = jac->kspupper;
1732: if (n) *n = nn;
1733: }
1734: return(0);
1735: }
1737: static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp)
1738: {
1739: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1743: if (!jac->schur) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitGetSubKSP()");
1744: PetscMalloc1(jac->nsplits,subksp);
1745: MatSchurComplementGetKSP(jac->schur,*subksp);
1747: (*subksp)[1] = jac->kspschur;
1748: if (n) *n = jac->nsplits;
1749: return(0);
1750: }
1752: static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp)
1753: {
1754: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1755: PetscErrorCode ierr;
1756: PetscInt cnt = 0;
1757: PC_FieldSplitLink ilink = jac->head;
1760: PetscMalloc1(jac->nsplits,subksp);
1761: while (ilink) {
1762: (*subksp)[cnt++] = ilink->ksp;
1763: ilink = ilink->next;
1764: }
1765: 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);
1766: if (n) *n = jac->nsplits;
1767: return(0);
1768: }
1770: /*@C
1771: PCFieldSplitRestrictIS - Restricts the fieldsplit ISs to be within a given IS.
1773: Input Parameters:
1774: + pc - the preconditioner context
1775: - is - the index set that defines the indices to which the fieldsplit is to be restricted
1777: Level: advanced
1779: @*/
1780: PetscErrorCode PCFieldSplitRestrictIS(PC pc,IS isy)
1781: {
1787: PetscTryMethod(pc,"PCFieldSplitRestrictIS_C",(PC,IS),(pc,isy));
1788: return(0);
1789: }
1792: static PetscErrorCode PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy)
1793: {
1794: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1795: PetscErrorCode ierr;
1796: PC_FieldSplitLink ilink = jac->head, next;
1797: PetscInt localsize,size,sizez,i;
1798: const PetscInt *ind, *indz;
1799: PetscInt *indc, *indcz;
1800: PetscBool flg;
1803: ISGetLocalSize(isy,&localsize);
1804: MPI_Scan(&localsize,&size,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)isy));
1805: size -= localsize;
1806: while(ilink) {
1807: IS isrl,isr;
1808: PC subpc;
1809: ISEmbed(ilink->is, isy, PETSC_TRUE, &isrl);
1810: ISGetLocalSize(isrl,&localsize);
1811: PetscMalloc1(localsize,&indc);
1812: ISGetIndices(isrl,&ind);
1813: PetscArraycpy(indc,ind,localsize);
1814: ISRestoreIndices(isrl,&ind);
1815: ISDestroy(&isrl);
1816: for (i=0; i<localsize; i++) *(indc+i) += size;
1817: ISCreateGeneral(PetscObjectComm((PetscObject)isy),localsize,indc,PETSC_OWN_POINTER,&isr);
1818: PetscObjectReference((PetscObject)isr);
1819: ISDestroy(&ilink->is);
1820: ilink->is = isr;
1821: PetscObjectReference((PetscObject)isr);
1822: ISDestroy(&ilink->is_col);
1823: ilink->is_col = isr;
1824: ISDestroy(&isr);
1825: KSPGetPC(ilink->ksp, &subpc);
1826: PetscObjectTypeCompare((PetscObject)subpc,PCFIELDSPLIT,&flg);
1827: if(flg) {
1828: IS iszl,isz;
1829: MPI_Comm comm;
1830: ISGetLocalSize(ilink->is,&localsize);
1831: comm = PetscObjectComm((PetscObject)ilink->is);
1832: ISEmbed(isy, ilink->is, PETSC_TRUE, &iszl);
1833: MPI_Scan(&localsize,&sizez,1,MPIU_INT,MPI_SUM,comm);
1834: sizez -= localsize;
1835: ISGetLocalSize(iszl,&localsize);
1836: PetscMalloc1(localsize,&indcz);
1837: ISGetIndices(iszl,&indz);
1838: PetscArraycpy(indcz,indz,localsize);
1839: ISRestoreIndices(iszl,&indz);
1840: ISDestroy(&iszl);
1841: for (i=0; i<localsize; i++) *(indcz+i) += sizez;
1842: ISCreateGeneral(comm,localsize,indcz,PETSC_OWN_POINTER,&isz);
1843: PCFieldSplitRestrictIS(subpc,isz);
1844: ISDestroy(&isz);
1845: }
1846: next = ilink->next;
1847: ilink = next;
1848: }
1849: jac->isrestrict = PETSC_TRUE;
1850: return(0);
1851: }
1853: static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is)
1854: {
1855: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1856: PetscErrorCode ierr;
1857: PC_FieldSplitLink ilink, next = jac->head;
1858: char prefix[128];
1861: if (jac->splitdefined) {
1862: PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1863: return(0);
1864: }
1865: PetscNew(&ilink);
1866: if (splitname) {
1867: PetscStrallocpy(splitname,&ilink->splitname);
1868: } else {
1869: PetscMalloc1(8,&ilink->splitname);
1870: PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);
1871: }
1872: 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 */
1873: PetscObjectReference((PetscObject)is);
1874: ISDestroy(&ilink->is);
1875: ilink->is = is;
1876: PetscObjectReference((PetscObject)is);
1877: ISDestroy(&ilink->is_col);
1878: ilink->is_col = is;
1879: ilink->next = NULL;
1880: KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1881: KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);
1882: PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1883: KSPSetType(ilink->ksp,KSPPREONLY);
1884: PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);
1886: PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);
1887: KSPSetOptionsPrefix(ilink->ksp,prefix);
1889: if (!next) {
1890: jac->head = ilink;
1891: ilink->previous = NULL;
1892: } else {
1893: while (next->next) {
1894: next = next->next;
1895: }
1896: next->next = ilink;
1897: ilink->previous = next;
1898: }
1899: jac->nsplits++;
1900: return(0);
1901: }
1903: /*@C
1904: PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner
1906: Logically Collective on PC
1908: Input Parameters:
1909: + pc - the preconditioner context
1910: . splitname - name of this split, if NULL the number of the split is used
1911: . n - the number of fields in this split
1912: - fields - the fields in this split
1914: Level: intermediate
1916: Notes:
1917: Use PCFieldSplitSetIS() to set a completely general set of indices as a field.
1919: The PCFieldSplitSetFields() is for defining fields as strided blocks. For example, if the block
1920: 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
1921: 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x....
1922: where the numbered entries indicate what is in the field.
1924: This function is called once per split (it creates a new split each time). Solve options
1925: for this split will be available under the prefix -fieldsplit_SPLITNAME_.
1927: Developer Note: This routine does not actually create the IS representing the split, that is delayed
1928: until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be
1929: available when this routine is called.
1931: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS()
1933: @*/
1934: PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1935: {
1941: if (n < 1) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname);
1943: PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt*,const PetscInt*),(pc,splitname,n,fields,fields_col));
1944: return(0);
1945: }
1947: /*@
1948: PCFieldSplitSetDiagUseAmat - set flag indicating whether to extract diagonal blocks from Amat (rather than Pmat)
1950: Logically Collective on PC
1952: Input Parameters:
1953: + pc - the preconditioner object
1954: - flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from
1956: Options Database:
1957: . -pc_fieldsplit_diag_use_amat
1959: Level: intermediate
1961: .seealso: PCFieldSplitGetDiagUseAmat(), PCFieldSplitSetOffDiagUseAmat(), PCFIELDSPLIT
1963: @*/
1964: PetscErrorCode PCFieldSplitSetDiagUseAmat(PC pc,PetscBool flg)
1965: {
1966: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1967: PetscBool isfs;
1972: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1973: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1974: jac->diag_use_amat = flg;
1975: return(0);
1976: }
1978: /*@
1979: PCFieldSplitGetDiagUseAmat - get the flag indicating whether to extract diagonal blocks from Amat (rather than Pmat)
1981: Logically Collective on PC
1983: Input Parameters:
1984: . pc - the preconditioner object
1986: Output Parameters:
1987: . flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from
1990: Level: intermediate
1992: .seealso: PCFieldSplitSetDiagUseAmat(), PCFieldSplitGetOffDiagUseAmat(), PCFIELDSPLIT
1994: @*/
1995: PetscErrorCode PCFieldSplitGetDiagUseAmat(PC pc,PetscBool *flg)
1996: {
1997: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1998: PetscBool isfs;
2004: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2005: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
2006: *flg = jac->diag_use_amat;
2007: return(0);
2008: }
2010: /*@
2011: PCFieldSplitSetOffDiagUseAmat - set flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat)
2013: Logically Collective on PC
2015: Input Parameters:
2016: + pc - the preconditioner object
2017: - flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from
2019: Options Database:
2020: . -pc_fieldsplit_off_diag_use_amat
2022: Level: intermediate
2024: .seealso: PCFieldSplitGetOffDiagUseAmat(), PCFieldSplitSetDiagUseAmat(), PCFIELDSPLIT
2026: @*/
2027: PetscErrorCode PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg)
2028: {
2029: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2030: PetscBool isfs;
2035: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2036: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
2037: jac->offdiag_use_amat = flg;
2038: return(0);
2039: }
2041: /*@
2042: PCFieldSplitGetOffDiagUseAmat - get the flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat)
2044: Logically Collective on PC
2046: Input Parameters:
2047: . pc - the preconditioner object
2049: Output Parameters:
2050: . flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from
2053: Level: intermediate
2055: .seealso: PCFieldSplitSetOffDiagUseAmat(), PCFieldSplitGetDiagUseAmat(), PCFIELDSPLIT
2057: @*/
2058: PetscErrorCode PCFieldSplitGetOffDiagUseAmat(PC pc,PetscBool *flg)
2059: {
2060: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2061: PetscBool isfs;
2067: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2068: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
2069: *flg = jac->offdiag_use_amat;
2070: return(0);
2071: }
2075: /*@C
2076: PCFieldSplitSetIS - Sets the exact elements for field
2078: Logically Collective on PC
2080: Input Parameters:
2081: + pc - the preconditioner context
2082: . splitname - name of this split, if NULL the number of the split is used
2083: - is - the index set that defines the vector elements in this field
2086: Notes:
2087: Use PCFieldSplitSetFields(), for fields defined by strided types.
2089: This function is called once per split (it creates a new split each time). Solve options
2090: for this split will be available under the prefix -fieldsplit_SPLITNAME_.
2092: Level: intermediate
2094: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize()
2096: @*/
2097: PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is)
2098: {
2105: PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));
2106: return(0);
2107: }
2109: /*@C
2110: PCFieldSplitGetIS - Retrieves the elements for a field as an IS
2112: Logically Collective on PC
2114: Input Parameters:
2115: + pc - the preconditioner context
2116: - splitname - name of this split
2118: Output Parameter:
2119: - is - the index set that defines the vector elements in this field, or NULL if the field is not found
2121: Level: intermediate
2123: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS()
2125: @*/
2126: PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is)
2127: {
2134: {
2135: PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;
2136: PC_FieldSplitLink ilink = jac->head;
2137: PetscBool found;
2139: *is = NULL;
2140: while (ilink) {
2141: PetscStrcmp(ilink->splitname, splitname, &found);
2142: if (found) {
2143: *is = ilink->is;
2144: break;
2145: }
2146: ilink = ilink->next;
2147: }
2148: }
2149: return(0);
2150: }
2152: /*@C
2153: PCFieldSplitGetISByIndex - Retrieves the elements for a given index field as an IS
2155: Logically Collective on PC
2157: Input Parameters:
2158: + pc - the preconditioner context
2159: - index - index of this split
2161: Output Parameter:
2162: - is - the index set that defines the vector elements in this field
2164: Level: intermediate
2166: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitGetIS(), PCFieldSplitSetIS()
2168: @*/
2169: PetscErrorCode PCFieldSplitGetISByIndex(PC pc,PetscInt index,IS *is)
2170: {
2174: if (index < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",index);
2177: {
2178: PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;
2179: PC_FieldSplitLink ilink = jac->head;
2180: PetscInt i = 0;
2181: if (index >= jac->nsplits) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",index,jac->nsplits);
2183: while (i < index) {
2184: ilink = ilink->next;
2185: ++i;
2186: }
2187: PCFieldSplitGetIS(pc,ilink->splitname,is);
2188: }
2189: return(0);
2190: }
2192: /*@
2193: PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the
2194: fieldsplit preconditioner. If not set the matrix block size is used.
2196: Logically Collective on PC
2198: Input Parameters:
2199: + pc - the preconditioner context
2200: - bs - the block size
2202: Level: intermediate
2204: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields()
2206: @*/
2207: PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs)
2208: {
2214: PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));
2215: return(0);
2216: }
2218: /*@C
2219: PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits
2221: Collective on KSP
2223: Input Parameter:
2224: . pc - the preconditioner context
2226: Output Parameters:
2227: + n - the number of splits
2228: - subksp - the array of KSP contexts
2230: Note:
2231: After PCFieldSplitGetSubKSP() the array of KSPs is to be freed by the user with PetscFree()
2232: (not the KSP just the array that contains them).
2234: You must call PCSetUp() before calling PCFieldSplitGetSubKSP().
2236: If the fieldsplit is of type PC_COMPOSITE_SCHUR, it returns the KSP object used inside the
2237: Schur complement and the KSP object used to iterate over the Schur complement.
2238: To access all the KSP objects used in PC_COMPOSITE_SCHUR, use PCFieldSplitSchurGetSubKSP().
2240: If the fieldsplit is of type PC_COMPOSITE_GKB, it returns the KSP object used to solve the
2241: inner linear system defined by the matrix H in each loop.
2243: Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs.
2244: You can call PCFieldSplitGetSubKSP(pc,n,PETSC_NULL_KSP,ierr) to determine how large the
2245: KSP array must be.
2248: Level: advanced
2250: .seealso: PCFIELDSPLIT
2251: @*/
2252: PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
2253: {
2259: PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));
2260: return(0);
2261: }
2263: /*@C
2264: PCFieldSplitSchurGetSubKSP - Gets the KSP contexts used inside the Schur complement based PCFIELDSPLIT
2266: Collective on KSP
2268: Input Parameter:
2269: . pc - the preconditioner context
2271: Output Parameters:
2272: + n - the number of splits
2273: - subksp - the array of KSP contexts
2275: Note:
2276: After PCFieldSplitSchurGetSubKSP() the array of KSPs is to be freed by the user with PetscFree()
2277: (not the KSP just the array that contains them).
2279: You must call PCSetUp() before calling PCFieldSplitSchurGetSubKSP().
2281: If the fieldsplit type is of type PC_COMPOSITE_SCHUR, it returns (in order)
2282: - the KSP used for the (1,1) block
2283: - the KSP used for the Schur complement (not the one used for the interior Schur solver)
2284: - the KSP used for the (1,1) block in the upper triangular factor (if different from that of the (1,1) block).
2286: It returns a null array if the fieldsplit is not of type PC_COMPOSITE_SCHUR; in this case, you should use PCFieldSplitGetSubKSP().
2288: Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs.
2289: You can call PCFieldSplitSchurGetSubKSP(pc,n,PETSC_NULL_KSP,ierr) to determine how large the
2290: KSP array must be.
2292: Level: advanced
2294: .seealso: PCFIELDSPLIT
2295: @*/
2296: PetscErrorCode PCFieldSplitSchurGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
2297: {
2303: PetscUseMethod(pc,"PCFieldSplitSchurGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));
2304: return(0);
2305: }
2307: /*@
2308: PCFieldSplitSetSchurPre - Indicates from what operator the preconditioner is constructucted for the Schur complement.
2309: The default is the A11 matrix.
2311: Collective on PC
2313: Input Parameters:
2314: + pc - the preconditioner context
2315: . 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
2316: PC_FIELDSPLIT_SCHUR_PRE_SELFP, and PC_FIELDSPLIT_SCHUR_PRE_FULL
2317: - userpre - matrix to use for preconditioning, or NULL
2319: Options Database:
2320: + -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11. See notes for meaning of various arguments
2321: - -fieldsplit_1_pc_type <pctype> - the preconditioner algorithm that is used to construct the preconditioner from the operator
2323: Notes:
2324: $ If ptype is
2325: $ a11 - the preconditioner for the Schur complement is generated from the block diagonal part of the preconditioner
2326: $ matrix associated with the Schur complement (i.e. A11), not the Schur complement matrix
2327: $ self - the preconditioner for the Schur complement is generated from the symbolic representation of the Schur complement matrix:
2328: $ The only preconditioner that currently works with this symbolic respresentation matrix object is the PCLSC
2329: $ preconditioner
2330: $ user - the preconditioner for the Schur complement is generated from the user provided matrix (pre argument
2331: $ to this function).
2332: $ selfp - the preconditioning for the Schur complement is generated from an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01
2333: $ This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be
2334: $ lumped before extracting the diagonal using the additional option -fieldsplit_1_mat_schur_complement_ainv_type lump
2335: $ full - the preconditioner for the Schur complement is generated from the exact Schur complement matrix representation computed internally by PCFIELDSPLIT (this is expensive)
2336: $ useful mostly as a test that the Schur complement approach can work for your problem
2338: When solving a saddle point problem, where the A11 block is identically zero, using a11 as the ptype only makes sense
2339: with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and
2340: -fieldsplit_1_pc_type lsc which uses the least squares commutator to compute a preconditioner for the Schur complement.
2342: Level: intermediate
2344: .seealso: PCFieldSplitGetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType,
2345: MatSchurComplementSetAinvType(), PCLSC
2347: @*/
2348: PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
2349: {
2354: PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));
2355: return(0);
2356: }
2358: PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) {return PCFieldSplitSetSchurPre(pc,ptype,pre);} /* Deprecated name */
2360: /*@
2361: PCFieldSplitGetSchurPre - For Schur complement fieldsplit, determine how the Schur complement will be
2362: preconditioned. See PCFieldSplitSetSchurPre() for details.
2364: Logically Collective on PC
2366: Input Parameters:
2367: . pc - the preconditioner context
2369: Output Parameters:
2370: + ptype - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11, PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_PRE_USER
2371: - userpre - matrix to use for preconditioning (with PC_FIELDSPLIT_PRE_USER), or NULL
2373: Level: intermediate
2375: .seealso: PCFieldSplitSetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC
2377: @*/
2378: PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
2379: {
2384: PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre));
2385: return(0);
2386: }
2388: /*@
2389: PCFieldSplitSchurGetS - extract the MatSchurComplement object used by this PC in case it needs to be configured separately
2391: Not collective
2393: Input Parameter:
2394: . pc - the preconditioner context
2396: Output Parameter:
2397: . S - the Schur complement matrix
2399: Notes:
2400: This matrix should not be destroyed using MatDestroy(); rather, use PCFieldSplitSchurRestoreS().
2402: Level: advanced
2404: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurRestoreS()
2406: @*/
2407: PetscErrorCode PCFieldSplitSchurGetS(PC pc,Mat *S)
2408: {
2410: const char* t;
2411: PetscBool isfs;
2412: PC_FieldSplit *jac;
2416: PetscObjectGetType((PetscObject)pc,&t);
2417: PetscStrcmp(t,PCFIELDSPLIT,&isfs);
2418: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
2419: jac = (PC_FieldSplit*)pc->data;
2420: if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
2421: if (S) *S = jac->schur;
2422: return(0);
2423: }
2425: /*@
2426: PCFieldSplitSchurRestoreS - restores the MatSchurComplement object used by this PC
2428: Not collective
2430: Input Parameters:
2431: + pc - the preconditioner context
2432: - S - the Schur complement matrix
2434: Level: advanced
2436: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurGetS()
2438: @*/
2439: PetscErrorCode PCFieldSplitSchurRestoreS(PC pc,Mat *S)
2440: {
2442: const char* t;
2443: PetscBool isfs;
2444: PC_FieldSplit *jac;
2448: PetscObjectGetType((PetscObject)pc,&t);
2449: PetscStrcmp(t,PCFIELDSPLIT,&isfs);
2450: if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
2451: jac = (PC_FieldSplit*)pc->data;
2452: if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
2453: if (!S || *S != jac->schur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten");
2454: return(0);
2455: }
2458: static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
2459: {
2460: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2464: jac->schurpre = ptype;
2465: if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) {
2466: MatDestroy(&jac->schur_user);
2467: jac->schur_user = pre;
2468: PetscObjectReference((PetscObject)jac->schur_user);
2469: }
2470: return(0);
2471: }
2473: static PetscErrorCode PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
2474: {
2475: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2478: *ptype = jac->schurpre;
2479: *pre = jac->schur_user;
2480: return(0);
2481: }
2483: /*@
2484: PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain in the preconditioner
2486: Collective on PC
2488: Input Parameters:
2489: + pc - the preconditioner context
2490: - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default
2492: Options Database:
2493: . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full
2496: Level: intermediate
2498: Notes:
2499: The FULL factorization is
2501: $ (A B) = (1 0) (A 0) (1 Ainv*B) = L D U
2502: $ (C E) (C*Ainv 1) (0 S) (0 1 )
2504: 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,
2505: 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().
2507: $ If A and S are solved exactly
2508: $ *) FULL factorization is a direct solver.
2509: $ *) 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.
2510: $ *) 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.
2512: If the iteration count is very low, consider using KSPFGMRES or KSPGCR which can use one less preconditioner
2513: Section 1.5 Writing Application Codes with PETSc in this case. Note that the preconditioned operator may be highly non-normal, so such fast convergence may not be observed in practice.
2515: For symmetric problems in which A is positive definite and S is negative definite, DIAG can be used with KSPMINRES.
2517: 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).
2519: References:
2520: + 1. - Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000).
2521: - 2. - Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001).
2523: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCFieldSplitSetSchurScale()
2524: @*/
2525: PetscErrorCode PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype)
2526: {
2531: PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));
2532: return(0);
2533: }
2535: static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype)
2536: {
2537: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2540: jac->schurfactorization = ftype;
2541: return(0);
2542: }
2544: /*@
2545: PCFieldSplitSetSchurScale - Controls the sign flip of S for PC_FIELDSPLIT_SCHUR_FACT_DIAG.
2547: Collective on PC
2549: Input Parameters:
2550: + pc - the preconditioner context
2551: - scale - scaling factor for the Schur complement
2553: Options Database:
2554: . -pc_fieldsplit_schur_scale - default is -1.0
2556: Level: intermediate
2558: .seealso: PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurFactType, PCFieldSplitSetSchurScale()
2559: @*/
2560: PetscErrorCode PCFieldSplitSetSchurScale(PC pc,PetscScalar scale)
2561: {
2567: PetscTryMethod(pc,"PCFieldSplitSetSchurScale_C",(PC,PetscScalar),(pc,scale));
2568: return(0);
2569: }
2571: static PetscErrorCode PCFieldSplitSetSchurScale_FieldSplit(PC pc,PetscScalar scale)
2572: {
2573: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2576: jac->schurscale = scale;
2577: return(0);
2578: }
2580: /*@C
2581: PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement
2583: Collective on KSP
2585: Input Parameter:
2586: . pc - the preconditioner context
2588: Output Parameters:
2589: + A00 - the (0,0) block
2590: . A01 - the (0,1) block
2591: . A10 - the (1,0) block
2592: - A11 - the (1,1) block
2594: Level: advanced
2596: .seealso: PCFIELDSPLIT
2597: @*/
2598: PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11)
2599: {
2600: PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;
2604: if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach.");
2605: if (A00) *A00 = jac->pmat[0];
2606: if (A01) *A01 = jac->B;
2607: if (A10) *A10 = jac->C;
2608: if (A11) *A11 = jac->pmat[1];
2609: return(0);
2610: }
2612: /*@
2613: PCFieldSplitSetGKBTol - Sets the solver tolerance for the generalized Golub-Kahan bidiagonalization preconditioner.
2615: Collective on PC
2617: Notes:
2618: The generalized GKB algorithm uses a lower bound estimate of the error in energy norm as stopping criterion.
2619: It stops once the lower bound estimate undershoots the required solver tolerance. Although the actual error might be bigger than
2620: this estimate, the stopping criterion is satisfactory in practical cases [A13].
2622: [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013.
2624: Input Parameters:
2625: + pc - the preconditioner context
2626: - tolerance - the solver tolerance
2628: Options Database:
2629: . -pc_fieldsplit_gkb_tol - default is 1e-5
2631: Level: intermediate
2633: .seealso: PCFIELDSPLIT, PCFieldSplitSetGKBDelay(), PCFieldSplitSetGKBNu(), PCFieldSplitSetGKBMaxit()
2634: @*/
2635: PetscErrorCode PCFieldSplitSetGKBTol(PC pc,PetscReal tolerance)
2636: {
2642: PetscTryMethod(pc,"PCFieldSplitSetGKBTol_C",(PC,PetscReal),(pc,tolerance));
2643: return(0);
2644: }
2646: static PetscErrorCode PCFieldSplitSetGKBTol_FieldSplit(PC pc,PetscReal tolerance)
2647: {
2648: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2651: jac->gkbtol = tolerance;
2652: return(0);
2653: }
2656: /*@
2657: PCFieldSplitSetGKBMaxit - Sets the maximum number of iterations for the generalized Golub-Kahan bidiagonalization
2658: preconditioner.
2660: Collective on PC
2662: Input Parameters:
2663: + pc - the preconditioner context
2664: - maxit - the maximum number of iterations
2666: Options Database:
2667: . -pc_fieldsplit_gkb_maxit - default is 100
2669: Level: intermediate
2671: .seealso: PCFIELDSPLIT, PCFieldSplitSetGKBDelay(), PCFieldSplitSetGKBTol(), PCFieldSplitSetGKBNu()
2672: @*/
2673: PetscErrorCode PCFieldSplitSetGKBMaxit(PC pc,PetscInt maxit)
2674: {
2680: PetscTryMethod(pc,"PCFieldSplitSetGKBMaxit_C",(PC,PetscInt),(pc,maxit));
2681: return(0);
2682: }
2684: static PetscErrorCode PCFieldSplitSetGKBMaxit_FieldSplit(PC pc,PetscInt maxit)
2685: {
2686: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2689: jac->gkbmaxit = maxit;
2690: return(0);
2691: }
2693: /*@
2694: PCFieldSplitSetGKBDelay - Sets the delay in the lower bound error estimate in the generalized Golub-Kahan bidiagonalization
2695: preconditioner.
2697: Collective on PC
2699: Notes:
2700: 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
2701: is expressed as a truncated sum. The error at iteration k can only be measured at iteration (k + delay), and thus the algorithm needs
2702: 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
2704: [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013.
2706: Input Parameters:
2707: + pc - the preconditioner context
2708: - delay - the delay window in the lower bound estimate
2710: Options Database:
2711: . -pc_fieldsplit_gkb_delay - default is 5
2713: Level: intermediate
2715: .seealso: PCFIELDSPLIT, PCFieldSplitSetGKBNu(), PCFieldSplitSetGKBTol(), PCFieldSplitSetGKBMaxit()
2716: @*/
2717: PetscErrorCode PCFieldSplitSetGKBDelay(PC pc,PetscInt delay)
2718: {
2724: PetscTryMethod(pc,"PCFieldSplitSetGKBDelay_C",(PC,PetscInt),(pc,delay));
2725: return(0);
2726: }
2728: static PetscErrorCode PCFieldSplitSetGKBDelay_FieldSplit(PC pc,PetscInt delay)
2729: {
2730: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2733: jac->gkbdelay = delay;
2734: return(0);
2735: }
2737: /*@
2738: 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.
2740: Collective on PC
2742: Notes:
2743: 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,
2744: 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
2745: necessary to find a good balance in between the convergence of the inner and outer loop.
2747: 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.
2749: [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013.
2751: Input Parameters:
2752: + pc - the preconditioner context
2753: - nu - the shift parameter
2755: Options Database:
2756: . -pc_fieldsplit_gkb_nu - default is 1
2758: Level: intermediate
2760: .seealso: PCFIELDSPLIT, PCFieldSplitSetGKBDelay(), PCFieldSplitSetGKBTol(), PCFieldSplitSetGKBMaxit()
2761: @*/
2762: PetscErrorCode PCFieldSplitSetGKBNu(PC pc,PetscReal nu)
2763: {
2769: PetscTryMethod(pc,"PCFieldSplitSetGKBNu_C",(PC,PetscReal),(pc,nu));
2770: return(0);
2771: }
2773: static PetscErrorCode PCFieldSplitSetGKBNu_FieldSplit(PC pc,PetscReal nu)
2774: {
2775: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2778: jac->gkbnu = nu;
2779: return(0);
2780: }
2783: static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type)
2784: {
2785: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2789: jac->type = type;
2791: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",0);
2792: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",0);
2793: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",0);
2794: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);
2795: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",0);
2796: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBTol_C",0);
2797: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBMaxit_C",0);
2798: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBNu_C",0);
2799: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBDelay_C",0);
2801: if (type == PC_COMPOSITE_SCHUR) {
2802: pc->ops->apply = PCApply_FieldSplit_Schur;
2803: pc->ops->view = PCView_FieldSplit_Schur;
2805: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);
2806: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit);
2807: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit);
2808: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);
2809: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",PCFieldSplitSetSchurScale_FieldSplit);
2810: } else if (type == PC_COMPOSITE_GKB){
2811: pc->ops->apply = PCApply_FieldSplit_GKB;
2812: pc->ops->view = PCView_FieldSplit_GKB;
2814: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2815: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBTol_C",PCFieldSplitSetGKBTol_FieldSplit);
2816: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBMaxit_C",PCFieldSplitSetGKBMaxit_FieldSplit);
2817: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBNu_C",PCFieldSplitSetGKBNu_FieldSplit);
2818: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBDelay_C",PCFieldSplitSetGKBDelay_FieldSplit);
2819: } else {
2820: pc->ops->apply = PCApply_FieldSplit;
2821: pc->ops->view = PCView_FieldSplit;
2823: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2824: }
2825: return(0);
2826: }
2828: static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs)
2829: {
2830: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2833: if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs);
2834: 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);
2835: jac->bs = bs;
2836: return(0);
2837: }
2839: /*@
2840: PCFieldSplitSetType - Sets the type of fieldsplit preconditioner.
2842: Collective on PC
2844: Input Parameter:
2845: + pc - the preconditioner context
2846: - type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR
2848: Options Database Key:
2849: . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type
2851: Level: Intermediate
2853: .seealso: PCCompositeSetType()
2855: @*/
2856: PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type)
2857: {
2862: PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));
2863: return(0);
2864: }
2866: /*@
2867: PCFieldSplitGetType - Gets the type of fieldsplit preconditioner.
2869: Not collective
2871: Input Parameter:
2872: . pc - the preconditioner context
2874: Output Parameter:
2875: . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR
2877: Level: Intermediate
2879: .seealso: PCCompositeSetType()
2880: @*/
2881: PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type)
2882: {
2883: PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;
2888: *type = jac->type;
2889: return(0);
2890: }
2892: /*@
2893: PCFieldSplitSetDMSplits - Flags whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible.
2895: Logically Collective
2897: Input Parameters:
2898: + pc - the preconditioner context
2899: - flg - boolean indicating whether to use field splits defined by the DM
2901: Options Database Key:
2902: . -pc_fieldsplit_dm_splits
2904: Level: Intermediate
2906: .seealso: PCFieldSplitGetDMSplits()
2908: @*/
2909: PetscErrorCode PCFieldSplitSetDMSplits(PC pc,PetscBool flg)
2910: {
2911: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2912: PetscBool isfs;
2918: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2919: if (isfs) {
2920: jac->dm_splits = flg;
2921: }
2922: return(0);
2923: }
2926: /*@
2927: PCFieldSplitGetDMSplits - Returns flag indicating whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible.
2929: Logically Collective
2931: Input Parameter:
2932: . pc - the preconditioner context
2934: Output Parameter:
2935: . flg - boolean indicating whether to use field splits defined by the DM
2937: Level: Intermediate
2939: .seealso: PCFieldSplitSetDMSplits()
2941: @*/
2942: PetscErrorCode PCFieldSplitGetDMSplits(PC pc,PetscBool* flg)
2943: {
2944: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2945: PetscBool isfs;
2951: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2952: if (isfs) {
2953: if(flg) *flg = jac->dm_splits;
2954: }
2955: return(0);
2956: }
2958: /*@
2959: PCFieldSplitGetDetectSaddlePoint - Returns flag indicating whether PCFieldSplit will attempt to automatically determine fields based on zero diagonal entries.
2961: Logically Collective
2963: Input Parameter:
2964: . pc - the preconditioner context
2966: Output Parameter:
2967: . flg - boolean indicating whether to detect fields or not
2969: Level: Intermediate
2971: .seealso: PCFIELDSPLIT, PCFieldSplitSetDetectSaddlePoint()
2973: @*/
2974: PetscErrorCode PCFieldSplitGetDetectSaddlePoint(PC pc,PetscBool *flg)
2975: {
2976: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2979: *flg = jac->detect;
2980: return(0);
2981: }
2983: /*@
2984: PCFieldSplitSetDetectSaddlePoint - Sets flag indicating whether PCFieldSplit will attempt to automatically determine fields based on zero diagonal entries.
2986: Logically Collective
2988: Notes:
2989: Also sets the split type to PC_COMPOSITE_SCHUR (see PCFieldSplitSetType()) and the Schur preconditioner type to PC_FIELDSPLIT_SCHUR_PRE_SELF (see PCFieldSplitSetSchurPre()).
2991: Input Parameter:
2992: . pc - the preconditioner context
2994: Output Parameter:
2995: . flg - boolean indicating whether to detect fields or not
2997: Options Database Key:
2998: . -pc_fieldsplit_detect_saddle_point
3000: Level: Intermediate
3002: .seealso: PCFIELDSPLIT, PCFieldSplitSetDetectSaddlePoint(), PCFieldSplitSetType(), PCFieldSplitSetSchurPre()
3004: @*/
3005: PetscErrorCode PCFieldSplitSetDetectSaddlePoint(PC pc,PetscBool flg)
3006: {
3007: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
3011: jac->detect = flg;
3012: if (jac->detect) {
3013: PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);
3014: PCFieldSplitSetSchurPre(pc,PC_FIELDSPLIT_SCHUR_PRE_SELF,NULL);
3015: }
3016: return(0);
3017: }
3019: /* -------------------------------------------------------------------------------------*/
3020: /*MC
3021: PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual
3022: fields or groups of fields. See the users manual section "Solving Block Matrices" for more details.
3024: To set options on the solvers for each block append -fieldsplit_ to all the PC
3025: options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1
3027: To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP()
3028: and set the options directly on the resulting KSP object
3030: Level: intermediate
3032: Options Database Keys:
3033: + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split
3034: . -pc_fieldsplit_default - automatically add any fields to additional splits that have not
3035: been supplied explicitly by -pc_fieldsplit_%d_fields
3036: . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields)
3037: . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur,gkb> - type of relaxation or factorization splitting
3038: . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11; see PCFieldSplitSetSchurPre()
3039: . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> - set factorization type when using -pc_fieldsplit_type schur; see PCFieldSplitSetSchurFactType()
3040: - -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero diagonal and uses Schur complement with no preconditioner as the solver
3042: Options prefixes for inner solvers when using the Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ .
3043: For all other solvers they are -fieldsplit_%d_ for the dth field; use -fieldsplit_ for all fields.
3044: The options prefix for the inner solver when using the Golub-Kahan biadiagonalization preconditioner is -fieldsplit_0_
3046: Notes:
3047: Use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS()
3048: to define a field by an arbitrary collection of entries.
3050: If no fields are set the default is used. The fields are defined by entries strided by bs,
3051: beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(),
3052: if this is not called the block size defaults to the blocksize of the second matrix passed
3053: to KSPSetOperators()/PCSetOperators().
3055: $ For the Schur complement preconditioner if J = ( A00 A01 )
3056: $ ( A10 A11 )
3057: $ the preconditioner using full factorization is
3058: $ ( I -ksp(A00) A01 ) ( inv(A00) 0 ) ( I 0 )
3059: $ ( 0 I ) ( 0 ksp(S) ) ( -A10 ksp(A00) I )
3060: where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. S is the Schur complement
3061: $ S = A11 - A10 ksp(A00) A01
3062: 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
3063: in providing the SECOND split or 1 if not give). For PCFieldSplitGetSubKSP() when field number is 0,
3064: it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default
3065: A11 is used to construct a preconditioner for S, use PCFieldSplitSetSchurPre() for all the possible ways to construct the preconditioner for S.
3067: The factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above,
3068: diag gives
3069: $ ( inv(A00) 0 )
3070: $ ( 0 -ksp(S) )
3071: 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
3072: can be turned off with PCFieldSplitSetSchurScale() or by command line -pc_fieldsplit_schur_scale 1.0. The lower factorization is the inverse of
3073: $ ( A00 0 )
3074: $ ( A10 S )
3075: where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of
3076: $ ( A00 A01 )
3077: $ ( 0 S )
3078: where again the inverses of A00 and S are applied using KSPs.
3080: If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS
3081: is used automatically for a second block.
3083: The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1.
3084: Generally it should be used with the AIJ format.
3086: The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see,
3087: for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT
3088: inside a smoother resulting in "Distributive Smoothers".
3090: There is a nice discussion of block preconditioners in
3092: [El08] A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier-Stokes equations
3093: Howard Elman, V.E. Howle, John Shadid, Robert Shuttleworth, Ray Tuminaro, Journal of Computational Physics 227 (2008) 1790--1808
3094: http://chess.cs.umd.edu/~elman/papers/tax.pdf
3096: The Constrained Pressure Preconditioner (CPR) can be implemented using PCCOMPOSITE with PCGALERKIN. CPR first solves an R A P subsystem, updates the
3097: residual on all variables (PCCompositeSetType(pc,PC_COMPOSITE_MULTIPLICATIVE)), and then applies a simple ILU like preconditioner on all the variables.
3099: The generalized Golub-Kahan bidiagonalization preconditioner (gkb) can be applied to symmetric 2x2 block matrices of the shape
3100: $ ( A00 A01 )
3101: $ ( A01' 0 )
3102: 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'.
3103: 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_.
3105: [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013.
3107: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCLSC,
3108: PCFieldSplitGetSubKSP(), PCFieldSplitSchurGetSubKSP(), PCFieldSplitSetFields(),
3109: PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre(), PCFieldSplitSetSchurFactType(),
3110: MatSchurComplementSetAinvType(), PCFieldSplitSetSchurScale(), PCFieldSplitSetDetectSaddlePoint()
3111: M*/
3113: PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc)
3114: {
3116: PC_FieldSplit *jac;
3119: PetscNewLog(pc,&jac);
3121: jac->bs = -1;
3122: jac->nsplits = 0;
3123: jac->type = PC_COMPOSITE_MULTIPLICATIVE;
3124: jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */
3125: jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL;
3126: jac->schurscale = -1.0;
3127: jac->dm_splits = PETSC_TRUE;
3128: jac->detect = PETSC_FALSE;
3129: jac->gkbtol = 1e-5;
3130: jac->gkbdelay = 5;
3131: jac->gkbnu = 1;
3132: jac->gkbmaxit = 100;
3133: jac->gkbmonitor = PETSC_FALSE;
3135: pc->data = (void*)jac;
3137: pc->ops->apply = PCApply_FieldSplit;
3138: pc->ops->applytranspose = PCApplyTranspose_FieldSplit;
3139: pc->ops->setup = PCSetUp_FieldSplit;
3140: pc->ops->reset = PCReset_FieldSplit;
3141: pc->ops->destroy = PCDestroy_FieldSplit;
3142: pc->ops->setfromoptions = PCSetFromOptions_FieldSplit;
3143: pc->ops->view = PCView_FieldSplit;
3144: pc->ops->applyrichardson = 0;
3146: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurGetSubKSP_C",PCFieldSplitSchurGetSubKSP_FieldSplit);
3147: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
3148: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);
3149: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);
3150: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);
3151: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);
3152: PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",PCFieldSplitRestrictIS_FieldSplit);
3153: return(0);
3154: }