1 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/
2 #include <petsc/private/kspimpl.h> /* This is needed to provide the appropriate PETSC_EXTERN for KSP_Solve_FS ....*/
3 #include <petscdm.h>
4
5 const char *const PCFieldSplitSchurPreTypes[] = {"SELF","SELFP","A11","USER","FULL","PCFieldSplitSchurPreType","PC_FIELDSPLIT_SCHUR_PRE_",NULL};
6 const char *const PCFieldSplitSchurFactTypes[] = {"DIAG","LOWER","UPPER","FULL","PCFieldSplitSchurFactType","PC_FIELDSPLIT_SCHUR_FACT_",NULL};
7
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;
9
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 };
22
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;
34
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 */
46
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 */
57
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;
66
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 */
73
74 /* This helper is so that setting a user-provided preconditioning matrix is orthogonal to choosing to use it. This way the
75 * application-provided FormJacobian can provide this matrix without interfering with the user's (command-line) choices. */
FieldSplitSchurPre(PC_FieldSplit * jac)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 }
88
89
90 #include <petscdraw.h>
PCView_FieldSplit(PC pc,PetscViewer viewer)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;
98
99 PetscFunctionBegin;
100 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
101 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr);
102 if (iascii) {
103 if (jac->bs > 0) {
104 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs);CHKERRQ(ierr);
105 } else {
106 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D\n",PCCompositeTypes[jac->type],jac->nsplits);CHKERRQ(ierr);
107 }
108 if (pc->useAmat) {
109 ierr = PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");CHKERRQ(ierr);
110 }
111 if (jac->diag_use_amat) {
112 ierr = PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for diagonal blocks\n");CHKERRQ(ierr);
113 }
114 if (jac->offdiag_use_amat) {
115 ierr = PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for off-diagonal blocks\n");CHKERRQ(ierr);
116 }
117 ierr = PetscViewerASCIIPrintf(viewer," Solver info for each split is in the following KSP objects:\n");CHKERRQ(ierr);
118 for (i=0; i<jac->nsplits; i++) {
119 if (ilink->fields) {
120 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);CHKERRQ(ierr);
121 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr);
122 for (j=0; j<ilink->nfields; j++) {
123 if (j > 0) {
124 ierr = PetscViewerASCIIPrintf(viewer,",");CHKERRQ(ierr);
125 }
126 ierr = PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);CHKERRQ(ierr);
127 }
128 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr);
129 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr);
130 } else {
131 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);CHKERRQ(ierr);
132 }
133 ierr = KSPView(ilink->ksp,viewer);CHKERRQ(ierr);
134 ilink = ilink->next;
135 }
136 }
137
138 if (isdraw) {
139 PetscDraw draw;
140 PetscReal x,y,w,wd;
141
142 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr);
143 ierr = PetscDrawGetCurrentPoint(draw,&x,&y);CHKERRQ(ierr);
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 ierr = PetscDrawPushCurrentPoint(draw,x,y);CHKERRQ(ierr);
149 ierr = KSPView(ilink->ksp,viewer);CHKERRQ(ierr);
150 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr);
151 x += wd;
152 ilink = ilink->next;
153 }
154 }
155 PetscFunctionReturn(0);
156 }
157
PCView_FieldSplit_Schur(PC pc,PetscViewer viewer)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;
166
167 PetscFunctionBegin;
168 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
169 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr);
170 if (iascii) {
171 if (jac->bs > 0) {
172 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, blocksize = %D, factorization %s\n",jac->bs,PCFieldSplitSchurFactTypes[jac->schurfactorization]);CHKERRQ(ierr);
173 } else {
174 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, factorization %s\n",PCFieldSplitSchurFactTypes[jac->schurfactorization]);CHKERRQ(ierr);
175 }
176 if (pc->useAmat) {
177 ierr = PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");CHKERRQ(ierr);
178 }
179 switch (jac->schurpre) {
180 case PC_FIELDSPLIT_SCHUR_PRE_SELF:
181 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from S itself\n");CHKERRQ(ierr);
182 break;
183 case PC_FIELDSPLIT_SCHUR_PRE_SELFP:
184 ierr = MatSchurComplementGetAinvType(jac->schur,&atype);CHKERRQ(ierr);
185 ierr = 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 "));CHKERRQ(ierr);break;
186 case PC_FIELDSPLIT_SCHUR_PRE_A11:
187 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n");CHKERRQ(ierr);
188 break;
189 case PC_FIELDSPLIT_SCHUR_PRE_FULL:
190 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from the exact Schur complement\n");CHKERRQ(ierr);
191 break;
192 case PC_FIELDSPLIT_SCHUR_PRE_USER:
193 if (jac->schur_user) {
194 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from user provided matrix\n");CHKERRQ(ierr);
195 } else {
196 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n");CHKERRQ(ierr);
197 }
198 break;
199 default:
200 SETERRQ1(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Invalid Schur preconditioning type: %d", jac->schurpre);
201 }
202 ierr = PetscViewerASCIIPrintf(viewer," Split info:\n");CHKERRQ(ierr);
203 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
204 for (i=0; i<jac->nsplits; i++) {
205 if (ilink->fields) {
206 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);CHKERRQ(ierr);
207 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr);
208 for (j=0; j<ilink->nfields; j++) {
209 if (j > 0) {
210 ierr = PetscViewerASCIIPrintf(viewer,",");CHKERRQ(ierr);
211 }
212 ierr = PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);CHKERRQ(ierr);
213 }
214 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr);
215 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr);
216 } else {
217 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);CHKERRQ(ierr);
218 }
219 ilink = ilink->next;
220 }
221 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for A00 block\n");CHKERRQ(ierr);
222 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
223 if (jac->head) {
224 ierr = KSPView(jac->head->ksp,viewer);CHKERRQ(ierr);
225 } else {ierr = PetscViewerASCIIPrintf(viewer," not yet available\n");CHKERRQ(ierr);}
226 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
227 if (jac->head && jac->kspupper != jac->head->ksp) {
228 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for upper A00 in upper triangular factor \n");CHKERRQ(ierr);
229 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
230 if (jac->kspupper) {ierr = KSPView(jac->kspupper,viewer);CHKERRQ(ierr);}
231 else {ierr = PetscViewerASCIIPrintf(viewer," not yet available\n");CHKERRQ(ierr);}
232 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
233 }
234 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for S = A11 - A10 inv(A00) A01 \n");CHKERRQ(ierr);
235 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
236 if (jac->kspschur) {
237 ierr = KSPView(jac->kspschur,viewer);CHKERRQ(ierr);
238 } else {
239 ierr = PetscViewerASCIIPrintf(viewer," not yet available\n");CHKERRQ(ierr);
240 }
241 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
242 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
243 } else if (isdraw && jac->head) {
244 PetscDraw draw;
245 PetscReal x,y,w,wd,h;
246 PetscInt cnt = 2;
247 char str[32];
248
249 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr);
250 ierr = PetscDrawGetCurrentPoint(draw,&x,&y);CHKERRQ(ierr);
251 if (jac->kspupper != jac->head->ksp) cnt++;
252 w = 2*PetscMin(1.0 - x,x);
253 wd = w/(cnt + 1);
254
255 ierr = PetscSNPrintf(str,32,"Schur fact. %s",PCFieldSplitSchurFactTypes[jac->schurfactorization]);CHKERRQ(ierr);
256 ierr = PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);CHKERRQ(ierr);
257 y -= h;
258 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_USER && !jac->schur_user) {
259 ierr = PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[PC_FIELDSPLIT_SCHUR_PRE_A11]);CHKERRQ(ierr);
260 } else {
261 ierr = PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[jac->schurpre]);CHKERRQ(ierr);
262 }
263 ierr = PetscDrawStringBoxed(draw,x+wd*(cnt-1)/2.0,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);CHKERRQ(ierr);
264 y -= h;
265 x = x - wd*(cnt-1)/2.0;
266
267 ierr = PetscDrawPushCurrentPoint(draw,x,y);CHKERRQ(ierr);
268 ierr = KSPView(jac->head->ksp,viewer);CHKERRQ(ierr);
269 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr);
270 if (jac->kspupper != jac->head->ksp) {
271 x += wd;
272 ierr = PetscDrawPushCurrentPoint(draw,x,y);CHKERRQ(ierr);
273 ierr = KSPView(jac->kspupper,viewer);CHKERRQ(ierr);
274 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr);
275 }
276 x += wd;
277 ierr = PetscDrawPushCurrentPoint(draw,x,y);CHKERRQ(ierr);
278 ierr = KSPView(jac->kspschur,viewer);CHKERRQ(ierr);
279 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr);
280 }
281 PetscFunctionReturn(0);
282 }
283
PCView_FieldSplit_GKB(PC pc,PetscViewer viewer)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;
291
292 PetscFunctionBegin;
293 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
294 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr);
295 if (iascii) {
296 if (jac->bs > 0) {
297 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs);CHKERRQ(ierr);
298 } else {
299 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D\n",PCCompositeTypes[jac->type],jac->nsplits);CHKERRQ(ierr);
300 }
301 if (pc->useAmat) {
302 ierr = PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");CHKERRQ(ierr);
303 }
304 if (jac->diag_use_amat) {
305 ierr = PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for diagonal blocks\n");CHKERRQ(ierr);
306 }
307 if (jac->offdiag_use_amat) {
308 ierr = PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for off-diagonal blocks\n");CHKERRQ(ierr);
309 }
310
311 ierr = PetscViewerASCIIPrintf(viewer," Stopping tolerance=%.1e, delay in error estimate=%D, maximum iterations=%D\n",jac->gkbtol,jac->gkbdelay,jac->gkbmaxit);CHKERRQ(ierr);
312 ierr = PetscViewerASCIIPrintf(viewer," Solver info for H = A00 + nu*A01*A01' matrix:\n");CHKERRQ(ierr);
313 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
314
315 if (ilink->fields) {
316 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",0);CHKERRQ(ierr);
317 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr);
318 for (j=0; j<ilink->nfields; j++) {
319 if (j > 0) {
320 ierr = PetscViewerASCIIPrintf(viewer,",");CHKERRQ(ierr);
321 }
322 ierr = PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);CHKERRQ(ierr);
323 }
324 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr);
325 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr);
326 } else {
327 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",0);CHKERRQ(ierr);
328 }
329 ierr = KSPView(ilink->ksp,viewer);CHKERRQ(ierr);
330
331 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
332 }
333
334 if (isdraw) {
335 PetscDraw draw;
336 PetscReal x,y,w,wd;
337
338 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr);
339 ierr = PetscDrawGetCurrentPoint(draw,&x,&y);CHKERRQ(ierr);
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 ierr = PetscDrawPushCurrentPoint(draw,x,y);CHKERRQ(ierr);
345 ierr = KSPView(ilink->ksp,viewer);CHKERRQ(ierr);
346 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr);
347 x += wd;
348 ilink = ilink->next;
349 }
350 }
351 PetscFunctionReturn(0);
352 }
353
354
355 /* Precondition: jac->bs is set to a meaningful value */
PCFieldSplitSetRuntimeSplits_Private(PC pc)356 static PetscErrorCode PCFieldSplitSetRuntimeSplits_Private(PC pc)
357 {
358 PetscErrorCode ierr;
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];
363
364 PetscFunctionBegin;
365 ierr = PetscMalloc1(jac->bs,&ifields);CHKERRQ(ierr);
366 ierr = PetscMalloc1(jac->bs,&ifields_col);CHKERRQ(ierr);
367 for (i=0,flg=PETSC_TRUE;; i++) {
368 ierr = PetscSNPrintf(splitname,sizeof(splitname),"%D",i);CHKERRQ(ierr);
369 ierr = PetscSNPrintf(optionname,sizeof(optionname),"-pc_fieldsplit_%D_fields",i);CHKERRQ(ierr);
370 ierr = PetscSNPrintf(optionname_col,sizeof(optionname_col),"-pc_fieldsplit_%D_fields_col",i);CHKERRQ(ierr);
371 nfields = jac->bs;
372 nfields_col = jac->bs;
373 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,optionname,ifields,&nfields,&flg);CHKERRQ(ierr);
374 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,optionname_col,ifields_col,&nfields_col,&flg_col);CHKERRQ(ierr);
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 ierr = PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields);CHKERRQ(ierr);
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 ierr = PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields_col);CHKERRQ(ierr);
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 ierr = PetscFree(ifields);CHKERRQ(ierr);
392 ierr = PetscFree(ifields_col);CHKERRQ(ierr);
393 PetscFunctionReturn(0);
394 }
395
PCFieldSplitSetDefaults(PC pc)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;
403
404 PetscFunctionBegin;
405 /*
406 Kinda messy, but at least this now uses DMCreateFieldDecomposition().
407 Should probably be rewritten.
408 */
409 if (!ilink) {
410 ierr = PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_coupling",&coupling,NULL);CHKERRQ(ierr);
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;
418
419 ierr = DMCreateFieldDecomposition(pc->dm, &numFields, &fieldNames, &fields, &dms);CHKERRQ(ierr);
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;
426
427 ierr = PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%D_fields", i);CHKERRQ(ierr);
428 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix, optionname, ifields, &nfields, &flg);CHKERRQ(ierr);
429 if (!flg) break;
430 if (numFields > 128) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot currently support %d > 128 fields", numFields);
431 ierr = DMCreateSubDM(pc->dm, nfields, ifields, &compField, &subdm[i]);CHKERRQ(ierr);
432 if (nfields == 1) {
433 ierr = PCFieldSplitSetIS(pc, fieldNames[ifields[0]], compField);CHKERRQ(ierr);
434 } else {
435 ierr = PetscSNPrintf(splitname, sizeof(splitname), "%D", i);CHKERRQ(ierr);
436 ierr = PCFieldSplitSetIS(pc, splitname, compField);CHKERRQ(ierr);
437 }
438 ierr = ISDestroy(&compField);CHKERRQ(ierr);
439 for (j = 0; j < nfields; ++j) {
440 f = ifields[j];
441 ierr = PetscFree(fieldNames[f]);CHKERRQ(ierr);
442 ierr = ISDestroy(&fields[f]);CHKERRQ(ierr);
443 }
444 }
445 if (i == 0) {
446 for (f = 0; f < numFields; ++f) {
447 ierr = PCFieldSplitSetIS(pc, fieldNames[f], fields[f]);CHKERRQ(ierr);
448 ierr = PetscFree(fieldNames[f]);CHKERRQ(ierr);
449 ierr = ISDestroy(&fields[f]);CHKERRQ(ierr);
450 }
451 } else {
452 for (j=0; j<numFields; j++) {
453 ierr = DMDestroy(dms+j);CHKERRQ(ierr);
454 }
455 ierr = PetscFree(dms);CHKERRQ(ierr);
456 ierr = PetscMalloc1(i, &dms);CHKERRQ(ierr);
457 for (j = 0; j < i; ++j) dms[j] = subdm[j];
458 }
459 ierr = PetscFree(fieldNames);CHKERRQ(ierr);
460 ierr = PetscFree(fields);CHKERRQ(ierr);
461 if (dms) {
462 ierr = PetscInfo(pc, "Setting up physics based fieldsplit preconditioner using the embedded DM\n");CHKERRQ(ierr);
463 for (ilink = jac->head, i = 0; ilink; ilink = ilink->next, ++i) {
464 const char *prefix;
465 ierr = PetscObjectGetOptionsPrefix((PetscObject)(ilink->ksp),&prefix);CHKERRQ(ierr);
466 ierr = PetscObjectSetOptionsPrefix((PetscObject)(dms[i]), prefix);CHKERRQ(ierr);
467 ierr = KSPSetDM(ilink->ksp, dms[i]);CHKERRQ(ierr);
468 ierr = KSPSetDMActive(ilink->ksp, PETSC_FALSE);CHKERRQ(ierr);
469 {
470 PetscErrorCode (*func)(KSP,Mat,Mat,void*);
471 void *ctx;
472
473 ierr = DMKSPGetComputeOperators(pc->dm, &func, &ctx);CHKERRQ(ierr);
474 ierr = DMKSPSetComputeOperators(dms[i], func, ctx);CHKERRQ(ierr);
475 }
476 ierr = PetscObjectIncrementTabLevel((PetscObject)dms[i],(PetscObject)ilink->ksp,0);CHKERRQ(ierr);
477 ierr = DMDestroy(&dms[i]);CHKERRQ(ierr);
478 }
479 ierr = PetscFree(dms);CHKERRQ(ierr);
480 }
481 } else {
482 if (jac->bs <= 0) {
483 if (pc->pmat) {
484 ierr = MatGetBlockSize(pc->pmat,&jac->bs);CHKERRQ(ierr);
485 } else jac->bs = 1;
486 }
487
488 if (jac->detect) {
489 IS zerodiags,rest;
490 PetscInt nmin,nmax;
491
492 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr);
493 if (jac->diag_use_amat) {
494 ierr = MatFindZeroDiagonals(pc->mat,&zerodiags);CHKERRQ(ierr);
495 } else {
496 ierr = MatFindZeroDiagonals(pc->pmat,&zerodiags);CHKERRQ(ierr);
497 }
498 ierr = ISComplement(zerodiags,nmin,nmax,&rest);CHKERRQ(ierr);
499 ierr = PCFieldSplitSetIS(pc,"0",rest);CHKERRQ(ierr);
500 ierr = PCFieldSplitSetIS(pc,"1",zerodiags);CHKERRQ(ierr);
501 ierr = ISDestroy(&zerodiags);CHKERRQ(ierr);
502 ierr = ISDestroy(&rest);CHKERRQ(ierr);
503 } else if (coupling) {
504 IS coupling,rest;
505 PetscInt nmin,nmax;
506
507 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr);
508 if (jac->offdiag_use_amat) {
509 ierr = MatFindOffBlockDiagonalEntries(pc->mat,&coupling);CHKERRQ(ierr);
510 } else {
511 ierr = MatFindOffBlockDiagonalEntries(pc->pmat,&coupling);CHKERRQ(ierr);
512 }
513 ierr = ISCreateStride(PetscObjectComm((PetscObject)pc->mat),nmax-nmin,nmin,1,&rest);CHKERRQ(ierr);
514 ierr = ISSetIdentity(rest);CHKERRQ(ierr);
515 ierr = PCFieldSplitSetIS(pc,"0",rest);CHKERRQ(ierr);
516 ierr = PCFieldSplitSetIS(pc,"1",coupling);CHKERRQ(ierr);
517 ierr = ISDestroy(&coupling);CHKERRQ(ierr);
518 ierr = ISDestroy(&rest);CHKERRQ(ierr);
519 } else {
520 ierr = PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&fieldsplit_default,NULL);CHKERRQ(ierr);
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 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr);
525 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);}
526 }
527 if ((fieldsplit_default || !jac->splitdefined) && !jac->isrestrict) {
528 Mat M = pc->pmat;
529 PetscBool isnest;
530
531 ierr = PetscInfo(pc,"Using default splitting of fields\n");CHKERRQ(ierr);
532 ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATNEST,&isnest);CHKERRQ(ierr);
533 if (!isnest) {
534 M = pc->mat;
535 ierr = PetscObjectTypeCompare((PetscObject)pc->mat,MATNEST,&isnest);CHKERRQ(ierr);
536 }
537 if (isnest) {
538 IS *fields;
539 PetscInt nf;
540
541 ierr = MatNestGetSize(M,&nf,NULL);CHKERRQ(ierr);
542 ierr = PetscMalloc1(nf,&fields);CHKERRQ(ierr);
543 ierr = MatNestGetISs(M,fields,NULL);CHKERRQ(ierr);
544 for (i=0;i<nf;i++) {
545 ierr = PCFieldSplitSetIS(pc,NULL,fields[i]);CHKERRQ(ierr);
546 }
547 ierr = PetscFree(fields);CHKERRQ(ierr);
548 } else {
549 for (i=0; i<jac->bs; i++) {
550 char splitname[8];
551 ierr = PetscSNPrintf(splitname,sizeof(splitname),"%D",i);CHKERRQ(ierr);
552 ierr = PCFieldSplitSetFields(pc,splitname,1,&i,&i);CHKERRQ(ierr);
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;
563
564 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr);
565 ierr = ISComplement(ilink->is,nmin,nmax,&is2);CHKERRQ(ierr);
566 ierr = PCFieldSplitSetIS(pc,"1",is2);CHKERRQ(ierr);
567 ierr = ISDestroy(&is2);CHKERRQ(ierr);
568 } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Must provide at least two sets of fields to PCFieldSplit()");
569 }
570
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 PetscFunctionReturn(0);
573 }
574
MatGolubKahanComputeExplicitOperator(Mat A,Mat B,Mat C,Mat * H,PetscReal gkbnu)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;
580
581 PetscFunctionBegin;
582 ierr = MatHermitianTranspose(C,MAT_INITIAL_MATRIX,&T);CHKERRQ(ierr); /* Test if augmented matrix is symmetric */
583 ierr = MatAXPY(T,-1.0,B,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
584 ierr = MatNorm(T,NORM_1,&nrmT);CHKERRQ(ierr);
585 ierr = MatNorm(B,NORM_1,&nrmB);CHKERRQ(ierr);
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 ierr = MatHermitianTranspose(B,MAT_INITIAL_MATRIX,&BT);CHKERRQ(ierr);
594 ierr = MatMatMult(B,BT,MAT_INITIAL_MATRIX,PETSC_DEFAULT,H);CHKERRQ(ierr); /* H = A01*A01' */
595 ierr = MatAYPX(*H,gkbnu,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); /* H = A00 + nu*A01*A01' */
596
597 ierr = MatDestroy(&BT);CHKERRQ(ierr);
598 ierr = MatDestroy(&T);CHKERRQ(ierr);
599 PetscFunctionReturn(0);
600 }
601
602 PETSC_EXTERN PetscErrorCode PetscOptionsFindPairPrefix_Private(PetscOptions,const char pre[], const char name[],const char *value[],PetscBool *flg);
603
PCSetUp_FieldSplit(PC pc)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;
611
612 PetscFunctionBegin;
613 pc->failedreason = PC_NOERROR;
614 ierr = PCFieldSplitSetDefaults(pc);CHKERRQ(ierr);
615 nsplit = jac->nsplits;
616 ilink = jac->head;
617
618 /* get the matrices for each split */
619 if (!jac->issetup) {
620 PetscInt rstart,rend,nslots,bs;
621
622 jac->issetup = PETSC_TRUE;
623
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 ierr = MatGetOwnershipRange(pc->pmat,&rstart,&rend);CHKERRQ(ierr);
630 nslots = (rend - rstart)/bs;
631 for (i=0; i<nsplit; i++) {
632 if (jac->defaultsplit) {
633 ierr = ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+i,nsplit,&ilink->is);CHKERRQ(ierr);
634 ierr = ISDuplicate(ilink->is,&ilink->is_col);CHKERRQ(ierr);
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 ierr = PetscMalloc1(ilink->nfields*nslots,&ii);CHKERRQ(ierr);
639 ierr = PetscMalloc1(ilink->nfields*nslots,&jj);CHKERRQ(ierr);
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 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,ii,PETSC_OWN_POINTER,&ilink->is);CHKERRQ(ierr);
647 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,jj,PETSC_OWN_POINTER,&ilink->is_col);CHKERRQ(ierr);
648 ierr = ISSetBlockSize(ilink->is, nfields);CHKERRQ(ierr);
649 ierr = ISSetBlockSize(ilink->is_col, nfields);CHKERRQ(ierr);
650 } else {
651 ierr = ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields[0],bs,&ilink->is);CHKERRQ(ierr);
652 ierr = ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields_col[0],bs,&ilink->is_col);CHKERRQ(ierr);
653 }
654 }
655 ierr = ISSorted(ilink->is,&sorted);CHKERRQ(ierr);
656 if (ilink->is_col) { ierr = ISSorted(ilink->is_col,&sorted_col);CHKERRQ(ierr); }
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 }
661
662 ilink = jac->head;
663 if (!jac->pmat) {
664 Vec xtmp;
665
666 ierr = MatCreateVecs(pc->pmat,&xtmp,NULL);CHKERRQ(ierr);
667 ierr = PetscMalloc1(nsplit,&jac->pmat);CHKERRQ(ierr);
668 ierr = PetscMalloc2(nsplit,&jac->x,nsplit,&jac->y);CHKERRQ(ierr);
669 for (i=0; i<nsplit; i++) {
670 MatNullSpace sp;
671
672 /* Check for preconditioning matrix attached to IS */
673 ierr = PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &jac->pmat[i]);CHKERRQ(ierr);
674 if (jac->pmat[i]) {
675 ierr = PetscObjectReference((PetscObject) jac->pmat[i]);CHKERRQ(ierr);
676 if (jac->type == PC_COMPOSITE_SCHUR) {
677 jac->schur_user = jac->pmat[i];
678
679 ierr = PetscObjectReference((PetscObject) jac->schur_user);CHKERRQ(ierr);
680 }
681 } else {
682 const char *prefix;
683 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->pmat[i]);CHKERRQ(ierr);
684 ierr = KSPGetOptionsPrefix(ilink->ksp,&prefix);CHKERRQ(ierr);
685 ierr = MatSetOptionsPrefix(jac->pmat[i],prefix);CHKERRQ(ierr);
686 ierr = MatViewFromOptions(jac->pmat[i],NULL,"-mat_view");CHKERRQ(ierr);
687 }
688 /* create work vectors for each split */
689 ierr = MatCreateVecs(jac->pmat[i],&jac->x[i],&jac->y[i]);CHKERRQ(ierr);
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 ierr = VecScatterCreate(xtmp,ilink->is,jac->x[i],NULL,&ilink->sctx);CHKERRQ(ierr);
693 ierr = PetscObjectQuery((PetscObject) ilink->is, "nearnullspace", (PetscObject*) &sp);CHKERRQ(ierr);
694 if (sp) {
695 ierr = MatSetNearNullSpace(jac->pmat[i], sp);CHKERRQ(ierr);
696 }
697 ilink = ilink->next;
698 }
699 ierr = VecDestroy(&xtmp);CHKERRQ(ierr);
700 } else {
701 MatReuse scall;
702 if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
703 for (i=0; i<nsplit; i++) {
704 ierr = MatDestroy(&jac->pmat[i]);CHKERRQ(ierr);
705 }
706 scall = MAT_INITIAL_MATRIX;
707 } else scall = MAT_REUSE_MATRIX;
708
709 for (i=0; i<nsplit; i++) {
710 Mat pmat;
711
712 /* Check for preconditioning matrix attached to IS */
713 ierr = PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &pmat);CHKERRQ(ierr);
714 if (!pmat) {
715 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is_col,scall,&jac->pmat[i]);CHKERRQ(ierr);
716 }
717 ilink = ilink->next;
718 }
719 }
720 if (jac->diag_use_amat) {
721 ilink = jac->head;
722 if (!jac->mat) {
723 ierr = PetscMalloc1(nsplit,&jac->mat);CHKERRQ(ierr);
724 for (i=0; i<nsplit; i++) {
725 ierr = MatCreateSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->mat[i]);CHKERRQ(ierr);
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 ierr = MatDestroy(&jac->mat[i]);CHKERRQ(ierr);
733 }
734 scall = MAT_INITIAL_MATRIX;
735 } else scall = MAT_REUSE_MATRIX;
736
737 for (i=0; i<nsplit; i++) {
738 ierr = MatCreateSubMatrix(pc->mat,ilink->is,ilink->is_col,scall,&jac->mat[i]);CHKERRQ(ierr);
739 ilink = ilink->next;
740 }
741 }
742 } else {
743 jac->mat = jac->pmat;
744 }
745
746 /* Check for null space attached to IS */
747 ilink = jac->head;
748 for (i=0; i<nsplit; i++) {
749 MatNullSpace sp;
750
751 ierr = PetscObjectQuery((PetscObject) ilink->is, "nullspace", (PetscObject*) &sp);CHKERRQ(ierr);
752 if (sp) {
753 ierr = MatSetNullSpace(jac->mat[i], sp);CHKERRQ(ierr);
754 }
755 ilink = ilink->next;
756 }
757
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 ierr = PetscCalloc1(nsplit,&jac->Afield);CHKERRQ(ierr);
766 if (jac->offdiag_use_amat) {
767 ierr = MatCreateSubMatrix(pc->mat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);CHKERRQ(ierr);
768 } else {
769 ierr = MatCreateSubMatrix(pc->pmat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);CHKERRQ(ierr);
770 }
771 } else {
772 MatReuse scall;
773
774 if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
775 ierr = MatDestroy(&jac->Afield[1]);CHKERRQ(ierr);
776 scall = MAT_INITIAL_MATRIX;
777 } else scall = MAT_REUSE_MATRIX;
778
779 if (jac->offdiag_use_amat) {
780 ierr = MatCreateSubMatrix(pc->mat,ilink->next->is,ilink->is,scall,&jac->Afield[1]);CHKERRQ(ierr);
781 } else {
782 ierr = MatCreateSubMatrix(pc->pmat,ilink->next->is,ilink->is,scall,&jac->Afield[1]);CHKERRQ(ierr);
783 }
784 }
785 } else {
786 if (!jac->Afield) {
787 ierr = PetscMalloc1(nsplit,&jac->Afield);CHKERRQ(ierr);
788 for (i=0; i<nsplit; i++) {
789 if (jac->offdiag_use_amat) {
790 ierr = MatCreateSubMatrix(pc->mat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);CHKERRQ(ierr);
791 } else {
792 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);CHKERRQ(ierr);
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 ierr = MatDestroy(&jac->Afield[i]);CHKERRQ(ierr);
801 }
802 scall = MAT_INITIAL_MATRIX;
803 } else scall = MAT_REUSE_MATRIX;
804
805 for (i=0; i<nsplit; i++) {
806 if (jac->offdiag_use_amat) {
807 ierr = MatCreateSubMatrix(pc->mat,ilink->is,NULL,scall,&jac->Afield[i]);CHKERRQ(ierr);
808 } else {
809 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,NULL,scall,&jac->Afield[i]);CHKERRQ(ierr);
810 }
811 ilink = ilink->next;
812 }
813 }
814 }
815 }
816
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;
823
824 if (nsplit != 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"To use Schur complement preconditioner you must have exactly 2 fields");
825
826 /* If pc->mat is SPD, don't scale by -1 the Schur complement */
827 if (jac->schurscale == (PetscScalar)-1.0) {
828 ierr = MatGetOption(pc->pmat,MAT_SPD,&isspd);CHKERRQ(ierr);
829 jac->schurscale = (isspd == PETSC_TRUE) ? 1.0 : -1.0;
830 }
831
832 /* When extracting off-diagonal submatrices, we take complements from this range */
833 ierr = MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);CHKERRQ(ierr);
834
835 if (jac->schur) {
836 KSP kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper;
837 MatReuse scall;
838
839 if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
840 scall = MAT_INITIAL_MATRIX;
841 ierr = MatDestroy(&jac->B);CHKERRQ(ierr);
842 ierr = MatDestroy(&jac->C);CHKERRQ(ierr);
843 } else scall = MAT_REUSE_MATRIX;
844
845 ierr = MatSchurComplementGetKSP(jac->schur, &kspInner);CHKERRQ(ierr);
846 ilink = jac->head;
847 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr);
848 if (jac->offdiag_use_amat) {
849 ierr = MatCreateSubMatrix(pc->mat,ilink->is,ccis,scall,&jac->B);CHKERRQ(ierr);
850 } else {
851 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,ccis,scall,&jac->B);CHKERRQ(ierr);
852 }
853 ierr = ISDestroy(&ccis);CHKERRQ(ierr);
854 ilink = ilink->next;
855 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr);
856 if (jac->offdiag_use_amat) {
857 ierr = MatCreateSubMatrix(pc->mat,ilink->is,ccis,scall,&jac->C);CHKERRQ(ierr);
858 } else {
859 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,ccis,scall,&jac->C);CHKERRQ(ierr);
860 }
861 ierr = ISDestroy(&ccis);CHKERRQ(ierr);
862 ierr = MatSchurComplementUpdateSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);CHKERRQ(ierr);
863 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
864 ierr = MatDestroy(&jac->schurp);CHKERRQ(ierr);
865 ierr = MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);CHKERRQ(ierr);
866 }
867 if (kspA != kspInner) {
868 ierr = KSPSetOperators(kspA,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr);
869 }
870 if (kspUpper != kspA) {
871 ierr = KSPSetOperators(kspUpper,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr);
872 }
873 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));CHKERRQ(ierr);
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;
881
882 /* extract the A01 and A10 matrices */
883 ilink = jac->head;
884 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr);
885 if (jac->offdiag_use_amat) {
886 ierr = MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr);
887 } else {
888 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr);
889 }
890 ierr = ISDestroy(&ccis);CHKERRQ(ierr);
891 ilink = ilink->next;
892 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr);
893 if (jac->offdiag_use_amat) {
894 ierr = MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr);
895 } else {
896 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr);
897 }
898 ierr = ISDestroy(&ccis);CHKERRQ(ierr);
899
900 /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */
901 ierr = MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur);CHKERRQ(ierr);
902 ierr = MatSetType(jac->schur,MATSCHURCOMPLEMENT);CHKERRQ(ierr);
903 ierr = MatSchurComplementSetSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);CHKERRQ(ierr);
904 ierr = PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr);
905 ierr = MatSetOptionsPrefix(jac->schur,schurmatprefix);CHKERRQ(ierr);
906 ierr = MatSchurComplementGetKSP(jac->schur,&kspt);CHKERRQ(ierr);
907 ierr = KSPSetOptionsPrefix(kspt,schurmatprefix);CHKERRQ(ierr);
908
909 /* Note: this is not true in general */
910 ierr = MatGetNullSpace(jac->mat[1], &sp);CHKERRQ(ierr);
911 if (sp) {
912 ierr = MatSetNullSpace(jac->schur, sp);CHKERRQ(ierr);
913 }
914
915 ierr = PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);CHKERRQ(ierr);
916 ierr = PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);CHKERRQ(ierr);
917 if (flg) {
918 DM dmInner;
919 KSP kspInner;
920 PC pcInner;
921
922 ierr = MatSchurComplementGetKSP(jac->schur, &kspInner);CHKERRQ(ierr);
923 ierr = KSPReset(kspInner);CHKERRQ(ierr);
924 ierr = KSPSetOperators(kspInner,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr);
925 ierr = PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr);
926 /* Indent this deeper to emphasize the "inner" nature of this solver. */
927 ierr = PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject) pc, 2);CHKERRQ(ierr);
928 ierr = PetscObjectIncrementTabLevel((PetscObject)kspInner->pc, (PetscObject) pc, 2);CHKERRQ(ierr);
929 ierr = KSPSetOptionsPrefix(kspInner, schurprefix);CHKERRQ(ierr);
930
931 /* Set DM for new solver */
932 ierr = KSPGetDM(jac->head->ksp, &dmInner);CHKERRQ(ierr);
933 ierr = KSPSetDM(kspInner, dmInner);CHKERRQ(ierr);
934 ierr = KSPSetDMActive(kspInner, PETSC_FALSE);CHKERRQ(ierr);
935
936 /* Defaults to PCKSP as preconditioner */
937 ierr = KSPGetPC(kspInner, &pcInner);CHKERRQ(ierr);
938 ierr = PCSetType(pcInner, PCKSP);CHKERRQ(ierr);
939 ierr = PCKSPSetKSP(pcInner, jac->head->ksp);CHKERRQ(ierr);
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 ierr = KSPSetType(jac->head->ksp,KSPGMRES);CHKERRQ(ierr);
948 ierr = MatSchurComplementSetKSP(jac->schur,jac->head->ksp);CHKERRQ(ierr);
949 }
950 ierr = KSPSetOperators(jac->head->ksp,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr);
951 ierr = KSPSetFromOptions(jac->head->ksp);CHKERRQ(ierr);
952 ierr = MatSetFromOptions(jac->schur);CHKERRQ(ierr);
953
954 ierr = PetscObjectTypeCompare((PetscObject)jac->schur, MATSCHURCOMPLEMENT, &flg);CHKERRQ(ierr);
955 if (flg) { /* Need to do this otherwise PCSetUp_KSP will overwrite the amat of jac->head->ksp */
956 KSP kspInner;
957 PC pcInner;
958
959 ierr = MatSchurComplementGetKSP(jac->schur, &kspInner);CHKERRQ(ierr);
960 ierr = KSPGetPC(kspInner, &pcInner);CHKERRQ(ierr);
961 ierr = PetscObjectTypeCompare((PetscObject)pcInner, PCKSP, &flg);CHKERRQ(ierr);
962 if (flg) {
963 KSP ksp;
964
965 ierr = PCKSPGetKSP(pcInner, &ksp);CHKERRQ(ierr);
966 if (ksp == jac->head->ksp) {
967 ierr = PCSetUseAmat(pcInner, PETSC_TRUE);CHKERRQ(ierr);
968 }
969 }
970 }
971 ierr = PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);CHKERRQ(ierr);
972 ierr = PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);CHKERRQ(ierr);
973 if (flg) {
974 DM dmInner;
975
976 ierr = PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr);
977 ierr = KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper);CHKERRQ(ierr);
978 ierr = KSPSetErrorIfNotConverged(jac->kspupper,pc->erroriffailure);CHKERRQ(ierr);
979 ierr = KSPSetOptionsPrefix(jac->kspupper, schurprefix);CHKERRQ(ierr);
980 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->kspupper, (PetscObject) pc, 1);CHKERRQ(ierr);
981 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->kspupper->pc, (PetscObject) pc, 1);CHKERRQ(ierr);
982 ierr = KSPGetDM(jac->head->ksp, &dmInner);CHKERRQ(ierr);
983 ierr = KSPSetDM(jac->kspupper, dmInner);CHKERRQ(ierr);
984 ierr = KSPSetDMActive(jac->kspupper, PETSC_FALSE);CHKERRQ(ierr);
985 ierr = KSPSetFromOptions(jac->kspupper);CHKERRQ(ierr);
986 ierr = KSPSetOperators(jac->kspupper,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr);
987 ierr = VecDuplicate(jac->head->x, &jac->head->z);CHKERRQ(ierr);
988 } else {
989 jac->kspupper = jac->head->ksp;
990 ierr = PetscObjectReference((PetscObject) jac->head->ksp);CHKERRQ(ierr);
991 }
992
993 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
994 ierr = MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);CHKERRQ(ierr);
995 }
996 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&jac->kspschur);CHKERRQ(ierr);
997 ierr = KSPSetErrorIfNotConverged(jac->kspschur,pc->erroriffailure);CHKERRQ(ierr);
998 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);CHKERRQ(ierr);
999 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);CHKERRQ(ierr);
1000 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) {
1001 PC pcschur;
1002 ierr = KSPGetPC(jac->kspschur,&pcschur);CHKERRQ(ierr);
1003 ierr = PCSetType(pcschur,PCNONE);CHKERRQ(ierr);
1004 /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */
1005 } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) {
1006 ierr = MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user);CHKERRQ(ierr);
1007 }
1008 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));CHKERRQ(ierr);
1009 ierr = KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix);CHKERRQ(ierr);
1010 ierr = KSPSetOptionsPrefix(jac->kspschur, Dprefix);CHKERRQ(ierr);
1011 /* propagate DM */
1012 {
1013 DM sdm;
1014 ierr = KSPGetDM(jac->head->next->ksp, &sdm);CHKERRQ(ierr);
1015 if (sdm) {
1016 ierr = KSPSetDM(jac->kspschur, sdm);CHKERRQ(ierr);
1017 ierr = KSPSetDMActive(jac->kspschur, PETSC_FALSE);CHKERRQ(ierr);
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 ierr = KSPSetFromOptions(jac->kspschur);CHKERRQ(ierr);
1023 }
1024 ierr = MatAssemblyBegin(jac->schur,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1025 ierr = MatAssemblyEnd(jac->schur,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1026
1027 /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */
1028 ierr = PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_L",ilink->splitname);CHKERRQ(ierr);
1029 ierr = PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);
1030 if (!LSC_L) {ierr = PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);}
1031 if (LSC_L) {ierr = PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L);CHKERRQ(ierr);}
1032 ierr = PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_Lp",ilink->splitname);CHKERRQ(ierr);
1033 ierr = PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);
1034 if (!LSC_L) {ierr = PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);}
1035 if (LSC_L) {ierr = PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L);CHKERRQ(ierr);}
1036 } else if (jac->type == PC_COMPOSITE_GKB) {
1037 IS ccis;
1038 PetscInt rstart,rend;
1039
1040 if (nsplit != 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"To use GKB preconditioner you must have exactly 2 fields");
1041
1042 ilink = jac->head;
1043
1044 /* When extracting off-diagonal submatrices, we take complements from this range */
1045 ierr = MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);CHKERRQ(ierr);
1046
1047 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr);
1048 if (jac->offdiag_use_amat) {
1049 ierr = MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr);
1050 } else {
1051 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr);
1052 }
1053 ierr = ISDestroy(&ccis);CHKERRQ(ierr);
1054 /* Create work vectors for GKB algorithm */
1055 ierr = VecDuplicate(ilink->x,&jac->u);CHKERRQ(ierr);
1056 ierr = VecDuplicate(ilink->x,&jac->Hu);CHKERRQ(ierr);
1057 ierr = VecDuplicate(ilink->x,&jac->w2);CHKERRQ(ierr);
1058 ilink = ilink->next;
1059 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr);
1060 if (jac->offdiag_use_amat) {
1061 ierr = MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr);
1062 } else {
1063 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr);
1064 }
1065 ierr = ISDestroy(&ccis);CHKERRQ(ierr);
1066 /* Create work vectors for GKB algorithm */
1067 ierr = VecDuplicate(ilink->x,&jac->v);CHKERRQ(ierr);
1068 ierr = VecDuplicate(ilink->x,&jac->d);CHKERRQ(ierr);
1069 ierr = VecDuplicate(ilink->x,&jac->w1);CHKERRQ(ierr);
1070 ierr = MatGolubKahanComputeExplicitOperator(jac->mat[0],jac->B,jac->C,&jac->H,jac->gkbnu);CHKERRQ(ierr);
1071 ierr = PetscCalloc1(jac->gkbdelay,&jac->vecz);CHKERRQ(ierr);
1072
1073 ilink = jac->head;
1074 ierr = KSPSetOperators(ilink->ksp,jac->H,jac->H);CHKERRQ(ierr);
1075 if (!jac->suboptionsset) {ierr = KSPSetFromOptions(ilink->ksp);CHKERRQ(ierr);}
1076 /* Create gkb_monitor context */
1077 if (jac->gkbmonitor) {
1078 PetscInt tablevel;
1079 ierr = PetscViewerCreate(PETSC_COMM_WORLD,&jac->gkbviewer);CHKERRQ(ierr);
1080 ierr = PetscViewerSetType(jac->gkbviewer,PETSCVIEWERASCII);CHKERRQ(ierr);
1081 ierr = PetscObjectGetTabLevel((PetscObject)ilink->ksp,&tablevel);CHKERRQ(ierr);
1082 ierr = PetscViewerASCIISetTab(jac->gkbviewer,tablevel);CHKERRQ(ierr);
1083 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)ilink->ksp,1);CHKERRQ(ierr);
1084 }
1085 } else {
1086 /* set up the individual splits' PCs */
1087 i = 0;
1088 ilink = jac->head;
1089 while (ilink) {
1090 ierr = KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i]);CHKERRQ(ierr);
1091 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
1092 if (!jac->suboptionsset) {ierr = KSPSetFromOptions(ilink->ksp);CHKERRQ(ierr);}
1093 i++;
1094 ilink = ilink->next;
1095 }
1096 }
1097
1098 jac->suboptionsset = PETSC_TRUE;
1099 PetscFunctionReturn(0);
1100 }
1101
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))
1111
PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y)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;
1118
1119 PetscFunctionBegin;
1120 switch (jac->schurfactorization) {
1121 case PC_FIELDSPLIT_SCHUR_FACT_DIAG:
1122 /* [A00 0; 0 -S], positive definite, suitable for MINRES */
1123 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1124 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1125 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1126 ierr = PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr);
1127 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr);
1128 ierr = KSPCheckSolve(kspA,pc,ilinkA->y);CHKERRQ(ierr);
1129 ierr = PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr);
1130 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1131 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1132 ierr = PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr);
1133 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr);
1134 ierr = KSPCheckSolve(jac->kspschur,pc,ilinkD->y);CHKERRQ(ierr);
1135 ierr = PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr);
1136 ierr = VecScale(ilinkD->y,jac->schurscale);CHKERRQ(ierr);
1137 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1138 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1139 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1140 break;
1141 case PC_FIELDSPLIT_SCHUR_FACT_LOWER:
1142 /* [A00 0; A10 S], suitable for left preconditioning */
1143 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1144 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1145 ierr = PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr);
1146 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr);
1147 ierr = KSPCheckSolve(kspA,pc,ilinkA->y);CHKERRQ(ierr);
1148 ierr = PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr);
1149 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr);
1150 ierr = VecScale(ilinkD->x,-1.);CHKERRQ(ierr);
1151 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1152 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1153 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1154 ierr = PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr);
1155 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr);
1156 ierr = KSPCheckSolve(jac->kspschur,pc,ilinkD->y);CHKERRQ(ierr);
1157 ierr = PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr);
1158 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1159 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1160 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1161 break;
1162 case PC_FIELDSPLIT_SCHUR_FACT_UPPER:
1163 /* [A00 A01; 0 S], suitable for right preconditioning */
1164 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1165 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1166 ierr = PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr);
1167 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr);
1168 ierr = KSPCheckSolve(jac->kspschur,pc,ilinkD->y);CHKERRQ(ierr);
1169 ierr = PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr);
1170 ierr = VecScale(ilinkA->x,-1.);CHKERRQ(ierr);
1171 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1172 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1173 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1174 ierr = PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr);
1175 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr);
1176 ierr = KSPCheckSolve(kspA,pc,ilinkA->y);CHKERRQ(ierr);
1177 ierr = PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr);
1178 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1179 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1180 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
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 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1185 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1186 ierr = PetscLogEventBegin(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr);
1187 ierr = KSPSolve(kspLower,ilinkA->x,ilinkA->y);CHKERRQ(ierr);
1188 ierr = KSPCheckSolve(kspLower,pc,ilinkA->y);CHKERRQ(ierr);
1189 ierr = PetscLogEventEnd(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr);
1190 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr);
1191 ierr = VecScale(ilinkD->x,-1.0);CHKERRQ(ierr);
1192 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1193 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1194
1195 ierr = PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr);
1196 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr);
1197 ierr = KSPCheckSolve(jac->kspschur,pc,ilinkD->y);CHKERRQ(ierr);
1198 ierr = PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr);
1199 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1200
1201 if (kspUpper == kspA) {
1202 ierr = MatMult(jac->B,ilinkD->y,ilinkA->y);CHKERRQ(ierr);
1203 ierr = VecAXPY(ilinkA->x,-1.0,ilinkA->y);CHKERRQ(ierr);
1204 ierr = PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr);
1205 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr);
1206 ierr = KSPCheckSolve(kspA,pc,ilinkA->y);CHKERRQ(ierr);
1207 ierr = PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr);
1208 } else {
1209 ierr = PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr);
1210 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr);
1211 ierr = KSPCheckSolve(kspA,pc,ilinkA->y);CHKERRQ(ierr);
1212 ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr);
1213 ierr = PetscLogEventBegin(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);CHKERRQ(ierr);
1214 ierr = KSPSolve(kspUpper,ilinkA->x,ilinkA->z);CHKERRQ(ierr);
1215 ierr = KSPCheckSolve(kspUpper,pc,ilinkA->z);CHKERRQ(ierr);
1216 ierr = PetscLogEventEnd(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);CHKERRQ(ierr);
1217 ierr = VecAXPY(ilinkA->y,-1.0,ilinkA->z);CHKERRQ(ierr);
1218 }
1219 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1220 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1221 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1222 }
1223 PetscFunctionReturn(0);
1224 }
1225
PCApply_FieldSplit(PC pc,Vec x,Vec y)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;
1232
1233 PetscFunctionBegin;
1234 if (jac->type == PC_COMPOSITE_ADDITIVE) {
1235 if (jac->defaultsplit) {
1236 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr);
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 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr);
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 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr);
1241 while (ilink) {
1242 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr);
1243 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr);
1244 ierr = KSPCheckSolve(ilink->ksp,pc,ilink->y);CHKERRQ(ierr);
1245 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr);
1246 ilink = ilink->next;
1247 }
1248 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr);
1249 } else {
1250 ierr = VecSet(y,0.0);CHKERRQ(ierr);
1251 while (ilink) {
1252 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr);
1253 ilink = ilink->next;
1254 }
1255 }
1256 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) {
1257 ierr = VecSet(y,0.0);CHKERRQ(ierr);
1258 /* solve on first block for first block variables */
1259 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1260 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1261 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr);
1262 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr);
1263 ierr = KSPCheckSolve(ilink->ksp,pc,ilink->y);CHKERRQ(ierr);
1264 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr);
1265 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1266 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1267
1268 /* compute the residual only onto second block variables using first block variables */
1269 ierr = MatMult(jac->Afield[1],ilink->y,ilink->next->x);CHKERRQ(ierr);
1270 ilink = ilink->next;
1271 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr);
1272 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1273 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1274
1275 /* solve on second block variables */
1276 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr);
1277 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr);
1278 ierr = KSPCheckSolve(ilink->ksp,pc,ilink->y);CHKERRQ(ierr);
1279 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr);
1280 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1281 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1282 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1283 if (!jac->w1) {
1284 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr);
1285 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr);
1286 }
1287 ierr = VecSet(y,0.0);CHKERRQ(ierr);
1288 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr);
1289 cnt = 1;
1290 while (ilink->next) {
1291 ilink = ilink->next;
1292 /* compute the residual only over the part of the vector needed */
1293 ierr = MatMult(jac->Afield[cnt++],y,ilink->x);CHKERRQ(ierr);
1294 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr);
1295 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1296 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1297 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr);
1298 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr);
1299 ierr = KSPCheckSolve(ilink->ksp,pc,ilink->y);CHKERRQ(ierr);
1300 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr);
1301 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1302 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
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 ierr = MatMult(jac->Afield[cnt--],y,ilink->x);CHKERRQ(ierr);
1310 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr);
1311 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1312 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1313 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr);
1314 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr);
1315 ierr = KSPCheckSolve(ilink->ksp,pc,ilink->y);CHKERRQ(ierr);
1316 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr);
1317 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1318 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1319 }
1320 }
1321 } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type);
1322 PetscFunctionReturn(0);
1323 }
1324
1325
PCApply_FieldSplit_GKB(PC pc,Vec x,Vec y)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;
1336
1337 PetscFunctionBegin;
1338 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1339 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1340 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1341 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1342
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;
1350
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 ierr = VecScale(ilinkD->x,jac->gkbnu);CHKERRQ(ierr);
1356 ierr = MatMultAdd(jac->B,ilinkD->x,ilinkA->x,ilinkA->x);CHKERRQ(ierr); /* 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 }
1362
1363 /* Transform rhs from [q,tilde{b}] to [0,b] */
1364 ierr = PetscLogEventBegin(ilinkA->event,ksp,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr);
1365 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr);
1366 ierr = KSPCheckSolve(ksp,pc,ilinkA->y);CHKERRQ(ierr);
1367 ierr = PetscLogEventEnd(ilinkA->event,ksp,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr);
1368 ierr = MatMultHermitianTranspose(jac->B,ilinkA->y,work1);CHKERRQ(ierr);
1369 ierr = VecAXPBY(work1,1.0/nu,-1.0,ilinkD->x);CHKERRQ(ierr); /* c = b - B'*x */
1370
1371 /* First step of algorithm */
1372 ierr = VecNorm(work1,NORM_2,&beta);CHKERRQ(ierr); /* beta = sqrt(nu*c'*c)*/
1373 KSPCheckDot(ksp,beta);
1374 beta = PetscSqrtScalar(nu)*beta;
1375 ierr = VecAXPBY(v,nu/beta,0.0,work1);CHKERRQ(ierr); /* v = nu/beta *c */
1376 ierr = MatMult(jac->B,v,work2);CHKERRQ(ierr); /* u = H^{-1}*B*v */
1377 ierr = PetscLogEventBegin(ilinkA->event,ksp,work2,u,NULL);CHKERRQ(ierr);
1378 ierr = KSPSolve(ksp,work2,u);CHKERRQ(ierr);
1379 ierr = KSPCheckSolve(ksp,pc,u);CHKERRQ(ierr);
1380 ierr = PetscLogEventEnd(ilinkA->event,ksp,work2,u,NULL);CHKERRQ(ierr);
1381 ierr = MatMult(jac->H,u,Hu);CHKERRQ(ierr); /* alpha = u'*H*u */
1382 ierr = VecDot(Hu,u,&alpha);CHKERRQ(ierr);
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 ierr = VecScale(u,1.0/alpha);CHKERRQ(ierr);
1387 ierr = VecAXPBY(d,1.0/alpha,0.0,v);CHKERRQ(ierr); /* v = nu/beta *c */
1388
1389 z = beta/alpha;
1390 vecz[1] = z;
1391
1392 /* Computation of first iterate x(1) and p(1) */
1393 ierr = VecAXPY(ilinkA->y,z,u);CHKERRQ(ierr);
1394 ierr = VecCopy(d,ilinkD->y);CHKERRQ(ierr);
1395 ierr = VecScale(ilinkD->y,-z);CHKERRQ(ierr);
1396
1397 iterGKB = 1; lowbnd = 2*jac->gkbtol;
1398 if (jac->gkbmonitor) {
1399 ierr = PetscViewerASCIIPrintf(jac->gkbviewer,"%3D GKB Lower bound estimate %14.12e\n",iterGKB,lowbnd);CHKERRQ(ierr);
1400 }
1401
1402 while (iterGKB < jac->gkbmaxit && lowbnd > jac->gkbtol) {
1403 iterGKB += 1;
1404 ierr = MatMultHermitianTranspose(jac->B,u,work1);CHKERRQ(ierr); /* v <- nu*(B'*u-alpha/nu*v) */
1405 ierr = VecAXPBY(v,nu,-alpha,work1);CHKERRQ(ierr);
1406 ierr = VecNorm(v,NORM_2,&beta);CHKERRQ(ierr); /* beta = sqrt(nu)*v'*v */
1407 beta = beta/PetscSqrtScalar(nu);
1408 ierr = VecScale(v,1.0/beta);CHKERRQ(ierr);
1409 ierr = MatMult(jac->B,v,work2);CHKERRQ(ierr); /* u <- H^{-1}*(B*v-beta*H*u) */
1410 ierr = MatMult(jac->H,u,Hu);CHKERRQ(ierr);
1411 ierr = VecAXPY(work2,-beta,Hu);CHKERRQ(ierr);
1412 ierr = PetscLogEventBegin(ilinkA->event,ksp,work2,u,NULL);CHKERRQ(ierr);
1413 ierr = KSPSolve(ksp,work2,u);CHKERRQ(ierr);
1414 ierr = KSPCheckSolve(ksp,pc,u);CHKERRQ(ierr);
1415 ierr = PetscLogEventEnd(ilinkA->event,ksp,work2,u,NULL);CHKERRQ(ierr);
1416 ierr = MatMult(jac->H,u,Hu);CHKERRQ(ierr); /* alpha = u'*H*u */
1417 ierr = VecDot(Hu,u,&alpha);CHKERRQ(ierr);
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 ierr = VecScale(u,1.0/alpha);CHKERRQ(ierr);
1422
1423 z = -beta/alpha*z; /* z <- beta/alpha*z */
1424 vecz[0] = z;
1425
1426 /* Computation of new iterate x(i+1) and p(i+1) */
1427 ierr = VecAXPBY(d,1.0/alpha,-beta/alpha,v);CHKERRQ(ierr); /* d = (v-beta*d)/alpha */
1428 ierr = VecAXPY(ilinkA->y,z,u);CHKERRQ(ierr); /* r = r + z*u */
1429 ierr = VecAXPY(ilinkD->y,-z,d);CHKERRQ(ierr); /* p = p - z*d */
1430 ierr = MatMult(jac->H,ilinkA->y,Hu);CHKERRQ(ierr); /* ||u||_H = u'*H*u */
1431 ierr = VecDot(Hu,ilinkA->y,&nrmz2);CHKERRQ(ierr);
1432
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 }
1441
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 ierr = PetscViewerASCIIPrintf(jac->gkbviewer,"%3D GKB Lower bound estimate %14.12e\n",iterGKB,lowbnd);CHKERRQ(ierr);
1447 }
1448 }
1449
1450 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1451 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1452 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1453 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1454
1455 PetscFunctionReturn(0);
1456 }
1457
1458
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))
1468
PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y)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;
1475
1476 PetscFunctionBegin;
1477 if (jac->type == PC_COMPOSITE_ADDITIVE) {
1478 if (jac->defaultsplit) {
1479 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr);
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 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr);
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 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr);
1484 while (ilink) {
1485 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr);
1486 ierr = KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr);
1487 ierr = KSPCheckSolve(ilink->ksp,pc,ilink->y);CHKERRQ(ierr);
1488 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr);
1489 ilink = ilink->next;
1490 }
1491 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr);
1492 } else {
1493 ierr = VecSet(y,0.0);CHKERRQ(ierr);
1494 while (ilink) {
1495 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr);
1496 ilink = ilink->next;
1497 }
1498 }
1499 } else {
1500 if (!jac->w1) {
1501 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr);
1502 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr);
1503 }
1504 ierr = VecSet(y,0.0);CHKERRQ(ierr);
1505 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1506 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr);
1507 while (ilink->next) {
1508 ilink = ilink->next;
1509 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr);
1510 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr);
1511 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr);
1512 }
1513 while (ilink->previous) {
1514 ilink = ilink->previous;
1515 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr);
1516 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr);
1517 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr);
1518 }
1519 } else {
1520 while (ilink->next) { /* get to last entry in linked list */
1521 ilink = ilink->next;
1522 }
1523 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr);
1524 while (ilink->previous) {
1525 ilink = ilink->previous;
1526 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr);
1527 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr);
1528 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr);
1529 }
1530 }
1531 }
1532 PetscFunctionReturn(0);
1533 }
1534
PCReset_FieldSplit(PC pc)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;
1540
1541 PetscFunctionBegin;
1542 while (ilink) {
1543 ierr = KSPDestroy(&ilink->ksp);CHKERRQ(ierr);
1544 ierr = VecDestroy(&ilink->x);CHKERRQ(ierr);
1545 ierr = VecDestroy(&ilink->y);CHKERRQ(ierr);
1546 ierr = VecDestroy(&ilink->z);CHKERRQ(ierr);
1547 ierr = VecScatterDestroy(&ilink->sctx);CHKERRQ(ierr);
1548 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr);
1549 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr);
1550 ierr = PetscFree(ilink->splitname);CHKERRQ(ierr);
1551 ierr = PetscFree(ilink->fields);CHKERRQ(ierr);
1552 ierr = PetscFree(ilink->fields_col);CHKERRQ(ierr);
1553 next = ilink->next;
1554 ierr = PetscFree(ilink);CHKERRQ(ierr);
1555 ilink = next;
1556 }
1557 jac->head = NULL;
1558 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr);
1559 if (jac->mat && jac->mat != jac->pmat) {
1560 ierr = MatDestroyMatrices(jac->nsplits,&jac->mat);CHKERRQ(ierr);
1561 } else if (jac->mat) {
1562 jac->mat = NULL;
1563 }
1564 if (jac->pmat) {ierr = MatDestroyMatrices(jac->nsplits,&jac->pmat);CHKERRQ(ierr);}
1565 if (jac->Afield) {ierr = MatDestroyMatrices(jac->nsplits,&jac->Afield);CHKERRQ(ierr);}
1566 jac->nsplits = 0;
1567 ierr = VecDestroy(&jac->w1);CHKERRQ(ierr);
1568 ierr = VecDestroy(&jac->w2);CHKERRQ(ierr);
1569 ierr = MatDestroy(&jac->schur);CHKERRQ(ierr);
1570 ierr = MatDestroy(&jac->schurp);CHKERRQ(ierr);
1571 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr);
1572 ierr = KSPDestroy(&jac->kspschur);CHKERRQ(ierr);
1573 ierr = KSPDestroy(&jac->kspupper);CHKERRQ(ierr);
1574 ierr = MatDestroy(&jac->B);CHKERRQ(ierr);
1575 ierr = MatDestroy(&jac->C);CHKERRQ(ierr);
1576 ierr = MatDestroy(&jac->H);CHKERRQ(ierr);
1577 ierr = VecDestroy(&jac->u);CHKERRQ(ierr);
1578 ierr = VecDestroy(&jac->v);CHKERRQ(ierr);
1579 ierr = VecDestroy(&jac->Hu);CHKERRQ(ierr);
1580 ierr = VecDestroy(&jac->d);CHKERRQ(ierr);
1581 ierr = PetscFree(jac->vecz);CHKERRQ(ierr);
1582 ierr = PetscViewerDestroy(&jac->gkbviewer);CHKERRQ(ierr);
1583 jac->isrestrict = PETSC_FALSE;
1584 PetscFunctionReturn(0);
1585 }
1586
PCDestroy_FieldSplit(PC pc)1587 static PetscErrorCode PCDestroy_FieldSplit(PC pc)
1588 {
1589 PetscErrorCode ierr;
1590
1591 PetscFunctionBegin;
1592 ierr = PCReset_FieldSplit(pc);CHKERRQ(ierr);
1593 ierr = PetscFree(pc->data);CHKERRQ(ierr);
1594 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurGetSubKSP_C",NULL);CHKERRQ(ierr);
1595 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL);CHKERRQ(ierr);
1596 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL);CHKERRQ(ierr);
1597 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL);CHKERRQ(ierr);
1598 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL);CHKERRQ(ierr);
1599 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL);CHKERRQ(ierr);
1600 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL);CHKERRQ(ierr);
1601 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL);CHKERRQ(ierr);
1602 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL);CHKERRQ(ierr);
1603 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",NULL);CHKERRQ(ierr);
1604 PetscFunctionReturn(0);
1605 }
1606
PCSetFromOptions_FieldSplit(PetscOptionItems * PetscOptionsObject,PC pc)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;
1614
1615 PetscFunctionBegin;
1616 ierr = PetscOptionsHead(PetscOptionsObject,"FieldSplit options");CHKERRQ(ierr);
1617 ierr = PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL);CHKERRQ(ierr);
1618 ierr = PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);CHKERRQ(ierr);
1619 if (flg) {
1620 ierr = PCFieldSplitSetBlockSize(pc,bs);CHKERRQ(ierr);
1621 }
1622 jac->diag_use_amat = pc->useAmat;
1623 ierr = 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);CHKERRQ(ierr);
1624 jac->offdiag_use_amat = pc->useAmat;
1625 ierr = 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);CHKERRQ(ierr);
1626 ierr = PetscOptionsBool("-pc_fieldsplit_detect_saddle_point","Form 2-way split by detecting zero diagonal entries", "PCFieldSplitSetDetectSaddlePoint",jac->detect,&jac->detect,NULL);CHKERRQ(ierr);
1627 ierr = PCFieldSplitSetDetectSaddlePoint(pc,jac->detect);CHKERRQ(ierr); /* Sets split type and Schur PC type */
1628 ierr = PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);CHKERRQ(ierr);
1629 if (flg) {
1630 ierr = PCFieldSplitSetType(pc,ctype);CHKERRQ(ierr);
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 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr);
1637 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);}
1638 }
1639 if (jac->type == PC_COMPOSITE_SCHUR) {
1640 ierr = PetscOptionsGetEnum(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);CHKERRQ(ierr);
1641 if (flg) {ierr = PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");CHKERRQ(ierr);}
1642 ierr = 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);CHKERRQ(ierr);
1643 ierr = PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSetSchurPre",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL);CHKERRQ(ierr);
1644 ierr = PetscOptionsScalar("-pc_fieldsplit_schur_scale","Scale Schur complement","PCFieldSplitSetSchurScale",jac->schurscale,&jac->schurscale,NULL);CHKERRQ(ierr);
1645 } else if (jac->type == PC_COMPOSITE_GKB) {
1646 ierr = PetscOptionsReal("-pc_fieldsplit_gkb_tol","The tolerance for the lower bound stopping criterion","PCFieldSplitGKBTol",jac->gkbtol,&jac->gkbtol,NULL);CHKERRQ(ierr);
1647 ierr = PetscOptionsInt("-pc_fieldsplit_gkb_delay","The delay value for lower bound criterion","PCFieldSplitGKBDelay",jac->gkbdelay,&jac->gkbdelay,NULL);CHKERRQ(ierr);
1648 ierr = PetscOptionsReal("-pc_fieldsplit_gkb_nu","Parameter in augmented Lagrangian approach","PCFieldSplitGKBNu",jac->gkbnu,&jac->gkbnu,NULL);CHKERRQ(ierr);
1649 if (jac->gkbnu < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nu cannot be less than 0: value %f",jac->gkbnu);
1650 ierr = PetscOptionsInt("-pc_fieldsplit_gkb_maxit","Maximum allowed number of iterations","PCFieldSplitGKBMaxit",jac->gkbmaxit,&jac->gkbmaxit,NULL);CHKERRQ(ierr);
1651 ierr = PetscOptionsBool("-pc_fieldsplit_gkb_monitor","Prints number of GKB iterations and error","PCFieldSplitGKB",jac->gkbmonitor,&jac->gkbmonitor,NULL);CHKERRQ(ierr);
1652 }
1653 ierr = PetscOptionsTail();CHKERRQ(ierr);
1654 PetscFunctionReturn(0);
1655 }
1656
1657 /*------------------------------------------------------------------------------------*/
1658
PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt * fields,const PetscInt * fields_col)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;
1666
1667 PetscFunctionBegin;
1668 if (jac->splitdefined) {
1669 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr);
1670 PetscFunctionReturn(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 ierr = PetscNew(&ilink);CHKERRQ(ierr);
1677 if (splitname) {
1678 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr);
1679 } else {
1680 ierr = PetscMalloc1(3,&ilink->splitname);CHKERRQ(ierr);
1681 ierr = PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);CHKERRQ(ierr);
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 ierr = PetscMalloc1(n,&ilink->fields);CHKERRQ(ierr);
1685 ierr = PetscArraycpy(ilink->fields,fields,n);CHKERRQ(ierr);
1686 ierr = PetscMalloc1(n,&ilink->fields_col);CHKERRQ(ierr);
1687 ierr = PetscArraycpy(ilink->fields_col,fields_col,n);CHKERRQ(ierr);
1688
1689 ilink->nfields = n;
1690 ilink->next = NULL;
1691 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);CHKERRQ(ierr);
1692 ierr = KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);CHKERRQ(ierr);
1693 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr);
1694 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr);
1695 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr);
1696
1697 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);CHKERRQ(ierr);
1698 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr);
1699
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 PetscFunctionReturn(0);
1712 }
1713
PCFieldSplitSchurGetSubKSP_FieldSplit(PC pc,PetscInt * n,KSP ** subksp)1714 static PetscErrorCode PCFieldSplitSchurGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp)
1715 {
1716 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1717 PetscErrorCode ierr;
1718
1719 PetscFunctionBegin;
1720 *subksp = NULL;
1721 if (n) *n = 0;
1722 if (jac->type == PC_COMPOSITE_SCHUR) {
1723 PetscInt nn;
1724
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 ierr = PetscMalloc1(nn,subksp);CHKERRQ(ierr);
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 PetscFunctionReturn(0);
1735 }
1736
PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt * n,KSP ** subksp)1737 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp)
1738 {
1739 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1740 PetscErrorCode ierr;
1741
1742 PetscFunctionBegin;
1743 if (!jac->schur) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitGetSubKSP()");
1744 ierr = PetscMalloc1(jac->nsplits,subksp);CHKERRQ(ierr);
1745 ierr = MatSchurComplementGetKSP(jac->schur,*subksp);CHKERRQ(ierr);
1746
1747 (*subksp)[1] = jac->kspschur;
1748 if (n) *n = jac->nsplits;
1749 PetscFunctionReturn(0);
1750 }
1751
PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt * n,KSP ** subksp)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;
1758
1759 PetscFunctionBegin;
1760 ierr = PetscMalloc1(jac->nsplits,subksp);CHKERRQ(ierr);
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 PetscFunctionReturn(0);
1768 }
1769
1770 /*@C
1771 PCFieldSplitRestrictIS - Restricts the fieldsplit ISs to be within a given IS.
1772
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
1776
1777 Level: advanced
1778
1779 @*/
PCFieldSplitRestrictIS(PC pc,IS isy)1780 PetscErrorCode PCFieldSplitRestrictIS(PC pc,IS isy)
1781 {
1782 PetscErrorCode ierr;
1783
1784 PetscFunctionBegin;
1785 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1786 PetscValidHeaderSpecific(isy,IS_CLASSID,2);
1787 ierr = PetscTryMethod(pc,"PCFieldSplitRestrictIS_C",(PC,IS),(pc,isy));CHKERRQ(ierr);
1788 PetscFunctionReturn(0);
1789 }
1790
1791
PCFieldSplitRestrictIS_FieldSplit(PC pc,IS isy)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;
1801
1802 PetscFunctionBegin;
1803 ierr = ISGetLocalSize(isy,&localsize);CHKERRQ(ierr);
1804 ierr = MPI_Scan(&localsize,&size,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)isy));CHKERRQ(ierr);
1805 size -= localsize;
1806 while (ilink) {
1807 IS isrl,isr;
1808 PC subpc;
1809 ierr = ISEmbed(ilink->is, isy, PETSC_TRUE, &isrl);CHKERRQ(ierr);
1810 ierr = ISGetLocalSize(isrl,&localsize);CHKERRQ(ierr);
1811 ierr = PetscMalloc1(localsize,&indc);CHKERRQ(ierr);
1812 ierr = ISGetIndices(isrl,&ind);CHKERRQ(ierr);
1813 ierr = PetscArraycpy(indc,ind,localsize);CHKERRQ(ierr);
1814 ierr = ISRestoreIndices(isrl,&ind);CHKERRQ(ierr);
1815 ierr = ISDestroy(&isrl);CHKERRQ(ierr);
1816 for (i=0; i<localsize; i++) *(indc+i) += size;
1817 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)isy),localsize,indc,PETSC_OWN_POINTER,&isr);CHKERRQ(ierr);
1818 ierr = PetscObjectReference((PetscObject)isr);CHKERRQ(ierr);
1819 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr);
1820 ilink->is = isr;
1821 ierr = PetscObjectReference((PetscObject)isr);CHKERRQ(ierr);
1822 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr);
1823 ilink->is_col = isr;
1824 ierr = ISDestroy(&isr);CHKERRQ(ierr);
1825 ierr = KSPGetPC(ilink->ksp, &subpc);CHKERRQ(ierr);
1826 ierr = PetscObjectTypeCompare((PetscObject)subpc,PCFIELDSPLIT,&flg);CHKERRQ(ierr);
1827 if (flg) {
1828 IS iszl,isz;
1829 MPI_Comm comm;
1830 ierr = ISGetLocalSize(ilink->is,&localsize);CHKERRQ(ierr);
1831 comm = PetscObjectComm((PetscObject)ilink->is);
1832 ierr = ISEmbed(isy, ilink->is, PETSC_TRUE, &iszl);CHKERRQ(ierr);
1833 ierr = MPI_Scan(&localsize,&sizez,1,MPIU_INT,MPI_SUM,comm);CHKERRQ(ierr);
1834 sizez -= localsize;
1835 ierr = ISGetLocalSize(iszl,&localsize);CHKERRQ(ierr);
1836 ierr = PetscMalloc1(localsize,&indcz);CHKERRQ(ierr);
1837 ierr = ISGetIndices(iszl,&indz);CHKERRQ(ierr);
1838 ierr = PetscArraycpy(indcz,indz,localsize);CHKERRQ(ierr);
1839 ierr = ISRestoreIndices(iszl,&indz);CHKERRQ(ierr);
1840 ierr = ISDestroy(&iszl);CHKERRQ(ierr);
1841 for (i=0; i<localsize; i++) *(indcz+i) += sizez;
1842 ierr = ISCreateGeneral(comm,localsize,indcz,PETSC_OWN_POINTER,&isz);CHKERRQ(ierr);
1843 ierr = PCFieldSplitRestrictIS(subpc,isz);CHKERRQ(ierr);
1844 ierr = ISDestroy(&isz);CHKERRQ(ierr);
1845 }
1846 next = ilink->next;
1847 ilink = next;
1848 }
1849 jac->isrestrict = PETSC_TRUE;
1850 PetscFunctionReturn(0);
1851 }
1852
PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is)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];
1859
1860 PetscFunctionBegin;
1861 if (jac->splitdefined) {
1862 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr);
1863 PetscFunctionReturn(0);
1864 }
1865 ierr = PetscNew(&ilink);CHKERRQ(ierr);
1866 if (splitname) {
1867 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr);
1868 } else {
1869 ierr = PetscMalloc1(8,&ilink->splitname);CHKERRQ(ierr);
1870 ierr = PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);CHKERRQ(ierr);
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 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr);
1874 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr);
1875 ilink->is = is;
1876 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr);
1877 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr);
1878 ilink->is_col = is;
1879 ilink->next = NULL;
1880 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);CHKERRQ(ierr);
1881 ierr = KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);CHKERRQ(ierr);
1882 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr);
1883 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr);
1884 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr);
1885
1886 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);CHKERRQ(ierr);
1887 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr);
1888
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 PetscFunctionReturn(0);
1901 }
1902
1903 /*@C
1904 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner
1905
1906 Logically Collective on PC
1907
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
1913
1914 Level: intermediate
1915
1916 Notes:
1917 Use PCFieldSplitSetIS() to set a completely general set of indices as a field.
1918
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.
1923
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_.
1926
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.
1930
1931 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS()
1932
1933 @*/
PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt * fields,const PetscInt * fields_col)1934 PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1935 {
1936 PetscErrorCode ierr;
1937
1938 PetscFunctionBegin;
1939 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1940 PetscValidCharPointer(splitname,2);
1941 if (n < 1) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname);
1942 PetscValidIntPointer(fields,3);
1943 ierr = PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt*,const PetscInt*),(pc,splitname,n,fields,fields_col));CHKERRQ(ierr);
1944 PetscFunctionReturn(0);
1945 }
1946
1947 /*@
1948 PCFieldSplitSetDiagUseAmat - set flag indicating whether to extract diagonal blocks from Amat (rather than Pmat)
1949
1950 Logically Collective on PC
1951
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
1955
1956 Options Database:
1957 . -pc_fieldsplit_diag_use_amat
1958
1959 Level: intermediate
1960
1961 .seealso: PCFieldSplitGetDiagUseAmat(), PCFieldSplitSetOffDiagUseAmat(), PCFIELDSPLIT
1962
1963 @*/
PCFieldSplitSetDiagUseAmat(PC pc,PetscBool flg)1964 PetscErrorCode PCFieldSplitSetDiagUseAmat(PC pc,PetscBool flg)
1965 {
1966 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1967 PetscBool isfs;
1968 PetscErrorCode ierr;
1969
1970 PetscFunctionBegin;
1971 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1972 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr);
1973 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1974 jac->diag_use_amat = flg;
1975 PetscFunctionReturn(0);
1976 }
1977
1978 /*@
1979 PCFieldSplitGetDiagUseAmat - get the flag indicating whether to extract diagonal blocks from Amat (rather than Pmat)
1980
1981 Logically Collective on PC
1982
1983 Input Parameters:
1984 . pc - the preconditioner object
1985
1986 Output Parameters:
1987 . flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from
1988
1989
1990 Level: intermediate
1991
1992 .seealso: PCFieldSplitSetDiagUseAmat(), PCFieldSplitGetOffDiagUseAmat(), PCFIELDSPLIT
1993
1994 @*/
PCFieldSplitGetDiagUseAmat(PC pc,PetscBool * flg)1995 PetscErrorCode PCFieldSplitGetDiagUseAmat(PC pc,PetscBool *flg)
1996 {
1997 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1998 PetscBool isfs;
1999 PetscErrorCode ierr;
2000
2001 PetscFunctionBegin;
2002 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
2003 PetscValidBoolPointer(flg,2);
2004 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr);
2005 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
2006 *flg = jac->diag_use_amat;
2007 PetscFunctionReturn(0);
2008 }
2009
2010 /*@
2011 PCFieldSplitSetOffDiagUseAmat - set flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat)
2012
2013 Logically Collective on PC
2014
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
2018
2019 Options Database:
2020 . -pc_fieldsplit_off_diag_use_amat
2021
2022 Level: intermediate
2023
2024 .seealso: PCFieldSplitGetOffDiagUseAmat(), PCFieldSplitSetDiagUseAmat(), PCFIELDSPLIT
2025
2026 @*/
PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg)2027 PetscErrorCode PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg)
2028 {
2029 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2030 PetscBool isfs;
2031 PetscErrorCode ierr;
2032
2033 PetscFunctionBegin;
2034 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
2035 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr);
2036 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
2037 jac->offdiag_use_amat = flg;
2038 PetscFunctionReturn(0);
2039 }
2040
2041 /*@
2042 PCFieldSplitGetOffDiagUseAmat - get the flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat)
2043
2044 Logically Collective on PC
2045
2046 Input Parameters:
2047 . pc - the preconditioner object
2048
2049 Output Parameters:
2050 . flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from
2051
2052
2053 Level: intermediate
2054
2055 .seealso: PCFieldSplitSetOffDiagUseAmat(), PCFieldSplitGetDiagUseAmat(), PCFIELDSPLIT
2056
2057 @*/
PCFieldSplitGetOffDiagUseAmat(PC pc,PetscBool * flg)2058 PetscErrorCode PCFieldSplitGetOffDiagUseAmat(PC pc,PetscBool *flg)
2059 {
2060 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2061 PetscBool isfs;
2062 PetscErrorCode ierr;
2063
2064 PetscFunctionBegin;
2065 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
2066 PetscValidBoolPointer(flg,2);
2067 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr);
2068 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
2069 *flg = jac->offdiag_use_amat;
2070 PetscFunctionReturn(0);
2071 }
2072
2073
2074
2075 /*@C
2076 PCFieldSplitSetIS - Sets the exact elements for field
2077
2078 Logically Collective on PC
2079
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
2084
2085
2086 Notes:
2087 Use PCFieldSplitSetFields(), for fields defined by strided types.
2088
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_.
2091
2092 Level: intermediate
2093
2094 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize()
2095
2096 @*/
PCFieldSplitSetIS(PC pc,const char splitname[],IS is)2097 PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is)
2098 {
2099 PetscErrorCode ierr;
2100
2101 PetscFunctionBegin;
2102 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
2103 if (splitname) PetscValidCharPointer(splitname,2);
2104 PetscValidHeaderSpecific(is,IS_CLASSID,3);
2105 ierr = PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));CHKERRQ(ierr);
2106 PetscFunctionReturn(0);
2107 }
2108
2109 /*@C
2110 PCFieldSplitGetIS - Retrieves the elements for a field as an IS
2111
2112 Logically Collective on PC
2113
2114 Input Parameters:
2115 + pc - the preconditioner context
2116 - splitname - name of this split
2117
2118 Output Parameter:
2119 - is - the index set that defines the vector elements in this field, or NULL if the field is not found
2120
2121 Level: intermediate
2122
2123 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS()
2124
2125 @*/
PCFieldSplitGetIS(PC pc,const char splitname[],IS * is)2126 PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is)
2127 {
2128 PetscErrorCode ierr;
2129
2130 PetscFunctionBegin;
2131 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
2132 PetscValidCharPointer(splitname,2);
2133 PetscValidPointer(is,3);
2134 {
2135 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;
2136 PC_FieldSplitLink ilink = jac->head;
2137 PetscBool found;
2138
2139 *is = NULL;
2140 while (ilink) {
2141 ierr = PetscStrcmp(ilink->splitname, splitname, &found);CHKERRQ(ierr);
2142 if (found) {
2143 *is = ilink->is;
2144 break;
2145 }
2146 ilink = ilink->next;
2147 }
2148 }
2149 PetscFunctionReturn(0);
2150 }
2151
2152 /*@C
2153 PCFieldSplitGetISByIndex - Retrieves the elements for a given index field as an IS
2154
2155 Logically Collective on PC
2156
2157 Input Parameters:
2158 + pc - the preconditioner context
2159 - index - index of this split
2160
2161 Output Parameter:
2162 - is - the index set that defines the vector elements in this field
2163
2164 Level: intermediate
2165
2166 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitGetIS(), PCFieldSplitSetIS()
2167
2168 @*/
PCFieldSplitGetISByIndex(PC pc,PetscInt index,IS * is)2169 PetscErrorCode PCFieldSplitGetISByIndex(PC pc,PetscInt index,IS *is)
2170 {
2171 PetscErrorCode ierr;
2172
2173 PetscFunctionBegin;
2174 if (index < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",index);
2175 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
2176 PetscValidPointer(is,3);
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);
2182
2183 while (i < index) {
2184 ilink = ilink->next;
2185 ++i;
2186 }
2187 ierr = PCFieldSplitGetIS(pc,ilink->splitname,is);CHKERRQ(ierr);
2188 }
2189 PetscFunctionReturn(0);
2190 }
2191
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.
2195
2196 Logically Collective on PC
2197
2198 Input Parameters:
2199 + pc - the preconditioner context
2200 - bs - the block size
2201
2202 Level: intermediate
2203
2204 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields()
2205
2206 @*/
PCFieldSplitSetBlockSize(PC pc,PetscInt bs)2207 PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs)
2208 {
2209 PetscErrorCode ierr;
2210
2211 PetscFunctionBegin;
2212 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
2213 PetscValidLogicalCollectiveInt(pc,bs,2);
2214 ierr = PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));CHKERRQ(ierr);
2215 PetscFunctionReturn(0);
2216 }
2217
2218 /*@C
2219 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits
2220
2221 Collective on KSP
2222
2223 Input Parameter:
2224 . pc - the preconditioner context
2225
2226 Output Parameters:
2227 + n - the number of splits
2228 - subksp - the array of KSP contexts
2229
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).
2233
2234 You must call PCSetUp() before calling PCFieldSplitGetSubKSP().
2235
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().
2239
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.
2242
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.
2246
2247
2248 Level: advanced
2249
2250 .seealso: PCFIELDSPLIT
2251 @*/
PCFieldSplitGetSubKSP(PC pc,PetscInt * n,KSP * subksp[])2252 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
2253 {
2254 PetscErrorCode ierr;
2255
2256 PetscFunctionBegin;
2257 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
2258 if (n) PetscValidIntPointer(n,2);
2259 ierr = PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));CHKERRQ(ierr);
2260 PetscFunctionReturn(0);
2261 }
2262
2263 /*@C
2264 PCFieldSplitSchurGetSubKSP - Gets the KSP contexts used inside the Schur complement based PCFIELDSPLIT
2265
2266 Collective on KSP
2267
2268 Input Parameter:
2269 . pc - the preconditioner context
2270
2271 Output Parameters:
2272 + n - the number of splits
2273 - subksp - the array of KSP contexts
2274
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).
2278
2279 You must call PCSetUp() before calling PCFieldSplitSchurGetSubKSP().
2280
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).
2285
2286 It returns a null array if the fieldsplit is not of type PC_COMPOSITE_SCHUR; in this case, you should use PCFieldSplitGetSubKSP().
2287
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.
2291
2292 Level: advanced
2293
2294 .seealso: PCFIELDSPLIT
2295 @*/
PCFieldSplitSchurGetSubKSP(PC pc,PetscInt * n,KSP * subksp[])2296 PetscErrorCode PCFieldSplitSchurGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
2297 {
2298 PetscErrorCode ierr;
2299
2300 PetscFunctionBegin;
2301 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
2302 if (n) PetscValidIntPointer(n,2);
2303 ierr = PetscUseMethod(pc,"PCFieldSplitSchurGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));CHKERRQ(ierr);
2304 PetscFunctionReturn(0);
2305 }
2306
2307 /*@
2308 PCFieldSplitSetSchurPre - Indicates from what operator the preconditioner is constructucted for the Schur complement.
2309 The default is the A11 matrix.
2310
2311 Collective on PC
2312
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
2318
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
2322
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
2337
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.
2341
2342 Level: intermediate
2343
2344 .seealso: PCFieldSplitGetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType,
2345 MatSchurComplementSetAinvType(), PCLSC
2346
2347 @*/
PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)2348 PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
2349 {
2350 PetscErrorCode ierr;
2351
2352 PetscFunctionBegin;
2353 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
2354 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));CHKERRQ(ierr);
2355 PetscFunctionReturn(0);
2356 }
2357
PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)2358 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) {return PCFieldSplitSetSchurPre(pc,ptype,pre);} /* Deprecated name */
2359
2360 /*@
2361 PCFieldSplitGetSchurPre - For Schur complement fieldsplit, determine how the Schur complement will be
2362 preconditioned. See PCFieldSplitSetSchurPre() for details.
2363
2364 Logically Collective on PC
2365
2366 Input Parameters:
2367 . pc - the preconditioner context
2368
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
2372
2373 Level: intermediate
2374
2375 .seealso: PCFieldSplitSetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC
2376
2377 @*/
PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType * ptype,Mat * pre)2378 PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
2379 {
2380 PetscErrorCode ierr;
2381
2382 PetscFunctionBegin;
2383 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
2384 ierr = PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre));CHKERRQ(ierr);
2385 PetscFunctionReturn(0);
2386 }
2387
2388 /*@
2389 PCFieldSplitSchurGetS - extract the MatSchurComplement object used by this PC in case it needs to be configured separately
2390
2391 Not collective
2392
2393 Input Parameter:
2394 . pc - the preconditioner context
2395
2396 Output Parameter:
2397 . S - the Schur complement matrix
2398
2399 Notes:
2400 This matrix should not be destroyed using MatDestroy(); rather, use PCFieldSplitSchurRestoreS().
2401
2402 Level: advanced
2403
2404 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurRestoreS()
2405
2406 @*/
PCFieldSplitSchurGetS(PC pc,Mat * S)2407 PetscErrorCode PCFieldSplitSchurGetS(PC pc,Mat *S)
2408 {
2409 PetscErrorCode ierr;
2410 const char* t;
2411 PetscBool isfs;
2412 PC_FieldSplit *jac;
2413
2414 PetscFunctionBegin;
2415 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
2416 ierr = PetscObjectGetType((PetscObject)pc,&t);CHKERRQ(ierr);
2417 ierr = PetscStrcmp(t,PCFIELDSPLIT,&isfs);CHKERRQ(ierr);
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 PetscFunctionReturn(0);
2423 }
2424
2425 /*@
2426 PCFieldSplitSchurRestoreS - restores the MatSchurComplement object used by this PC
2427
2428 Not collective
2429
2430 Input Parameters:
2431 + pc - the preconditioner context
2432 - S - the Schur complement matrix
2433
2434 Level: advanced
2435
2436 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurGetS()
2437
2438 @*/
PCFieldSplitSchurRestoreS(PC pc,Mat * S)2439 PetscErrorCode PCFieldSplitSchurRestoreS(PC pc,Mat *S)
2440 {
2441 PetscErrorCode ierr;
2442 const char* t;
2443 PetscBool isfs;
2444 PC_FieldSplit *jac;
2445
2446 PetscFunctionBegin;
2447 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
2448 ierr = PetscObjectGetType((PetscObject)pc,&t);CHKERRQ(ierr);
2449 ierr = PetscStrcmp(t,PCFIELDSPLIT,&isfs);CHKERRQ(ierr);
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 PetscFunctionReturn(0);
2455 }
2456
2457
PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)2458 static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
2459 {
2460 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2461 PetscErrorCode ierr;
2462
2463 PetscFunctionBegin;
2464 jac->schurpre = ptype;
2465 if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) {
2466 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr);
2467 jac->schur_user = pre;
2468 ierr = PetscObjectReference((PetscObject)jac->schur_user);CHKERRQ(ierr);
2469 }
2470 PetscFunctionReturn(0);
2471 }
2472
PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType * ptype,Mat * pre)2473 static PetscErrorCode PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
2474 {
2475 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2476
2477 PetscFunctionBegin;
2478 *ptype = jac->schurpre;
2479 *pre = jac->schur_user;
2480 PetscFunctionReturn(0);
2481 }
2482
2483 /*@
2484 PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain in the preconditioner
2485
2486 Collective on PC
2487
2488 Input Parameters:
2489 + pc - the preconditioner context
2490 - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default
2491
2492 Options Database:
2493 . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full
2494
2495
2496 Level: intermediate
2497
2498 Notes:
2499 The FULL factorization is
2500
2501 $ (A B) = (1 0) (A 0) (1 Ainv*B) = L D U
2502 $ (C E) (C*Ainv 1) (0 S) (0 1)
2503
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().
2506
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.
2511
2512 If the iteration count is very low, consider using KSPFGMRES or KSPGCR which can use one less preconditioner
2513 application in this case. Note that the preconditioned operator may be highly non-normal, so such fast convergence may not be observed in practice.
2514
2515 For symmetric problems in which A is positive definite and S is negative definite, DIAG can be used with KSPMINRES.
2516
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).
2518
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).
2522
2523 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCFieldSplitSetSchurScale()
2524 @*/
PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype)2525 PetscErrorCode PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype)
2526 {
2527 PetscErrorCode ierr;
2528
2529 PetscFunctionBegin;
2530 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
2531 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));CHKERRQ(ierr);
2532 PetscFunctionReturn(0);
2533 }
2534
PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype)2535 static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype)
2536 {
2537 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2538
2539 PetscFunctionBegin;
2540 jac->schurfactorization = ftype;
2541 PetscFunctionReturn(0);
2542 }
2543
2544 /*@
2545 PCFieldSplitSetSchurScale - Controls the sign flip of S for PC_FIELDSPLIT_SCHUR_FACT_DIAG.
2546
2547 Collective on PC
2548
2549 Input Parameters:
2550 + pc - the preconditioner context
2551 - scale - scaling factor for the Schur complement
2552
2553 Options Database:
2554 . -pc_fieldsplit_schur_scale - default is -1.0
2555
2556 Level: intermediate
2557
2558 .seealso: PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurFactType, PCFieldSplitSetSchurScale()
2559 @*/
PCFieldSplitSetSchurScale(PC pc,PetscScalar scale)2560 PetscErrorCode PCFieldSplitSetSchurScale(PC pc,PetscScalar scale)
2561 {
2562 PetscErrorCode ierr;
2563
2564 PetscFunctionBegin;
2565 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
2566 PetscValidLogicalCollectiveScalar(pc,scale,2);
2567 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurScale_C",(PC,PetscScalar),(pc,scale));CHKERRQ(ierr);
2568 PetscFunctionReturn(0);
2569 }
2570
PCFieldSplitSetSchurScale_FieldSplit(PC pc,PetscScalar scale)2571 static PetscErrorCode PCFieldSplitSetSchurScale_FieldSplit(PC pc,PetscScalar scale)
2572 {
2573 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2574
2575 PetscFunctionBegin;
2576 jac->schurscale = scale;
2577 PetscFunctionReturn(0);
2578 }
2579
2580 /*@C
2581 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement
2582
2583 Collective on KSP
2584
2585 Input Parameter:
2586 . pc - the preconditioner context
2587
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
2593
2594 Level: advanced
2595
2596 .seealso: PCFIELDSPLIT
2597 @*/
PCFieldSplitGetSchurBlocks(PC pc,Mat * A00,Mat * A01,Mat * A10,Mat * A11)2598 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11)
2599 {
2600 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;
2601
2602 PetscFunctionBegin;
2603 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
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 PetscFunctionReturn(0);
2610 }
2611
2612 /*@
2613 PCFieldSplitSetGKBTol - Sets the solver tolerance for the generalized Golub-Kahan bidiagonalization preconditioner.
2614
2615 Collective on PC
2616
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].
2621
2622 [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013.
2623
2624 Input Parameters:
2625 + pc - the preconditioner context
2626 - tolerance - the solver tolerance
2627
2628 Options Database:
2629 . -pc_fieldsplit_gkb_tol - default is 1e-5
2630
2631 Level: intermediate
2632
2633 .seealso: PCFIELDSPLIT, PCFieldSplitSetGKBDelay(), PCFieldSplitSetGKBNu(), PCFieldSplitSetGKBMaxit()
2634 @*/
PCFieldSplitSetGKBTol(PC pc,PetscReal tolerance)2635 PetscErrorCode PCFieldSplitSetGKBTol(PC pc,PetscReal tolerance)
2636 {
2637 PetscErrorCode ierr;
2638
2639 PetscFunctionBegin;
2640 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
2641 PetscValidLogicalCollectiveReal(pc,tolerance,2);
2642 ierr = PetscTryMethod(pc,"PCFieldSplitSetGKBTol_C",(PC,PetscReal),(pc,tolerance));CHKERRQ(ierr);
2643 PetscFunctionReturn(0);
2644 }
2645
PCFieldSplitSetGKBTol_FieldSplit(PC pc,PetscReal tolerance)2646 static PetscErrorCode PCFieldSplitSetGKBTol_FieldSplit(PC pc,PetscReal tolerance)
2647 {
2648 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2649
2650 PetscFunctionBegin;
2651 jac->gkbtol = tolerance;
2652 PetscFunctionReturn(0);
2653 }
2654
2655
2656 /*@
2657 PCFieldSplitSetGKBMaxit - Sets the maximum number of iterations for the generalized Golub-Kahan bidiagonalization
2658 preconditioner.
2659
2660 Collective on PC
2661
2662 Input Parameters:
2663 + pc - the preconditioner context
2664 - maxit - the maximum number of iterations
2665
2666 Options Database:
2667 . -pc_fieldsplit_gkb_maxit - default is 100
2668
2669 Level: intermediate
2670
2671 .seealso: PCFIELDSPLIT, PCFieldSplitSetGKBDelay(), PCFieldSplitSetGKBTol(), PCFieldSplitSetGKBNu()
2672 @*/
PCFieldSplitSetGKBMaxit(PC pc,PetscInt maxit)2673 PetscErrorCode PCFieldSplitSetGKBMaxit(PC pc,PetscInt maxit)
2674 {
2675 PetscErrorCode ierr;
2676
2677 PetscFunctionBegin;
2678 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
2679 PetscValidLogicalCollectiveInt(pc,maxit,2);
2680 ierr = PetscTryMethod(pc,"PCFieldSplitSetGKBMaxit_C",(PC,PetscInt),(pc,maxit));CHKERRQ(ierr);
2681 PetscFunctionReturn(0);
2682 }
2683
PCFieldSplitSetGKBMaxit_FieldSplit(PC pc,PetscInt maxit)2684 static PetscErrorCode PCFieldSplitSetGKBMaxit_FieldSplit(PC pc,PetscInt maxit)
2685 {
2686 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2687
2688 PetscFunctionBegin;
2689 jac->gkbmaxit = maxit;
2690 PetscFunctionReturn(0);
2691 }
2692
2693 /*@
2694 PCFieldSplitSetGKBDelay - Sets the delay in the lower bound error estimate in the generalized Golub-Kahan bidiagonalization
2695 preconditioner.
2696
2697 Collective on PC
2698
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
2703
2704 [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013.
2705
2706 Input Parameters:
2707 + pc - the preconditioner context
2708 - delay - the delay window in the lower bound estimate
2709
2710 Options Database:
2711 . -pc_fieldsplit_gkb_delay - default is 5
2712
2713 Level: intermediate
2714
2715 .seealso: PCFIELDSPLIT, PCFieldSplitSetGKBNu(), PCFieldSplitSetGKBTol(), PCFieldSplitSetGKBMaxit()
2716 @*/
PCFieldSplitSetGKBDelay(PC pc,PetscInt delay)2717 PetscErrorCode PCFieldSplitSetGKBDelay(PC pc,PetscInt delay)
2718 {
2719 PetscErrorCode ierr;
2720
2721 PetscFunctionBegin;
2722 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
2723 PetscValidLogicalCollectiveInt(pc,delay,2);
2724 ierr = PetscTryMethod(pc,"PCFieldSplitSetGKBDelay_C",(PC,PetscInt),(pc,delay));CHKERRQ(ierr);
2725 PetscFunctionReturn(0);
2726 }
2727
PCFieldSplitSetGKBDelay_FieldSplit(PC pc,PetscInt delay)2728 static PetscErrorCode PCFieldSplitSetGKBDelay_FieldSplit(PC pc,PetscInt delay)
2729 {
2730 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2731
2732 PetscFunctionBegin;
2733 jac->gkbdelay = delay;
2734 PetscFunctionReturn(0);
2735 }
2736
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.
2739
2740 Collective on PC
2741
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.
2746
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.
2748
2749 [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013.
2750
2751 Input Parameters:
2752 + pc - the preconditioner context
2753 - nu - the shift parameter
2754
2755 Options Database:
2756 . -pc_fieldsplit_gkb_nu - default is 1
2757
2758 Level: intermediate
2759
2760 .seealso: PCFIELDSPLIT, PCFieldSplitSetGKBDelay(), PCFieldSplitSetGKBTol(), PCFieldSplitSetGKBMaxit()
2761 @*/
PCFieldSplitSetGKBNu(PC pc,PetscReal nu)2762 PetscErrorCode PCFieldSplitSetGKBNu(PC pc,PetscReal nu)
2763 {
2764 PetscErrorCode ierr;
2765
2766 PetscFunctionBegin;
2767 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
2768 PetscValidLogicalCollectiveReal(pc,nu,2);
2769 ierr = PetscTryMethod(pc,"PCFieldSplitSetGKBNu_C",(PC,PetscReal),(pc,nu));CHKERRQ(ierr);
2770 PetscFunctionReturn(0);
2771 }
2772
PCFieldSplitSetGKBNu_FieldSplit(PC pc,PetscReal nu)2773 static PetscErrorCode PCFieldSplitSetGKBNu_FieldSplit(PC pc,PetscReal nu)
2774 {
2775 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2776
2777 PetscFunctionBegin;
2778 jac->gkbnu = nu;
2779 PetscFunctionReturn(0);
2780 }
2781
2782
PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type)2783 static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type)
2784 {
2785 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2786 PetscErrorCode ierr;
2787
2788 PetscFunctionBegin;
2789 jac->type = type;
2790
2791 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",0);CHKERRQ(ierr);
2792 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",0);CHKERRQ(ierr);
2793 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",0);CHKERRQ(ierr);
2794 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);CHKERRQ(ierr);
2795 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",0);CHKERRQ(ierr);
2796 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBTol_C",0);CHKERRQ(ierr);
2797 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBMaxit_C",0);CHKERRQ(ierr);
2798 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBNu_C",0);CHKERRQ(ierr);
2799 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBDelay_C",0);CHKERRQ(ierr);
2800
2801 if (type == PC_COMPOSITE_SCHUR) {
2802 pc->ops->apply = PCApply_FieldSplit_Schur;
2803 pc->ops->view = PCView_FieldSplit_Schur;
2804
2805 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);CHKERRQ(ierr);
2806 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit);CHKERRQ(ierr);
2807 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit);CHKERRQ(ierr);
2808 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);CHKERRQ(ierr);
2809 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",PCFieldSplitSetSchurScale_FieldSplit);CHKERRQ(ierr);
2810 } else if (type == PC_COMPOSITE_GKB){
2811 pc->ops->apply = PCApply_FieldSplit_GKB;
2812 pc->ops->view = PCView_FieldSplit_GKB;
2813
2814 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr);
2815 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBTol_C",PCFieldSplitSetGKBTol_FieldSplit);CHKERRQ(ierr);
2816 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBMaxit_C",PCFieldSplitSetGKBMaxit_FieldSplit);CHKERRQ(ierr);
2817 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBNu_C",PCFieldSplitSetGKBNu_FieldSplit);CHKERRQ(ierr);
2818 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBDelay_C",PCFieldSplitSetGKBDelay_FieldSplit);CHKERRQ(ierr);
2819 } else {
2820 pc->ops->apply = PCApply_FieldSplit;
2821 pc->ops->view = PCView_FieldSplit;
2822
2823 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr);
2824 }
2825 PetscFunctionReturn(0);
2826 }
2827
PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs)2828 static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs)
2829 {
2830 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2831
2832 PetscFunctionBegin;
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 PetscFunctionReturn(0);
2837 }
2838
2839 /*@
2840 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner.
2841
2842 Collective on PC
2843
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
2847
2848 Options Database Key:
2849 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type
2850
2851 Level: Intermediate
2852
2853 .seealso: PCCompositeSetType()
2854
2855 @*/
PCFieldSplitSetType(PC pc,PCCompositeType type)2856 PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type)
2857 {
2858 PetscErrorCode ierr;
2859
2860 PetscFunctionBegin;
2861 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
2862 ierr = PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));CHKERRQ(ierr);
2863 PetscFunctionReturn(0);
2864 }
2865
2866 /*@
2867 PCFieldSplitGetType - Gets the type of fieldsplit preconditioner.
2868
2869 Not collective
2870
2871 Input Parameter:
2872 . pc - the preconditioner context
2873
2874 Output Parameter:
2875 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR
2876
2877 Level: Intermediate
2878
2879 .seealso: PCCompositeSetType()
2880 @*/
PCFieldSplitGetType(PC pc,PCCompositeType * type)2881 PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type)
2882 {
2883 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;
2884
2885 PetscFunctionBegin;
2886 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
2887 PetscValidIntPointer(type,2);
2888 *type = jac->type;
2889 PetscFunctionReturn(0);
2890 }
2891
2892 /*@
2893 PCFieldSplitSetDMSplits - Flags whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible.
2894
2895 Logically Collective
2896
2897 Input Parameters:
2898 + pc - the preconditioner context
2899 - flg - boolean indicating whether to use field splits defined by the DM
2900
2901 Options Database Key:
2902 . -pc_fieldsplit_dm_splits
2903
2904 Level: Intermediate
2905
2906 .seealso: PCFieldSplitGetDMSplits()
2907
2908 @*/
PCFieldSplitSetDMSplits(PC pc,PetscBool flg)2909 PetscErrorCode PCFieldSplitSetDMSplits(PC pc,PetscBool flg)
2910 {
2911 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2912 PetscBool isfs;
2913 PetscErrorCode ierr;
2914
2915 PetscFunctionBegin;
2916 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
2917 PetscValidLogicalCollectiveBool(pc,flg,2);
2918 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr);
2919 if (isfs) {
2920 jac->dm_splits = flg;
2921 }
2922 PetscFunctionReturn(0);
2923 }
2924
2925
2926 /*@
2927 PCFieldSplitGetDMSplits - Returns flag indicating whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible.
2928
2929 Logically Collective
2930
2931 Input Parameter:
2932 . pc - the preconditioner context
2933
2934 Output Parameter:
2935 . flg - boolean indicating whether to use field splits defined by the DM
2936
2937 Level: Intermediate
2938
2939 .seealso: PCFieldSplitSetDMSplits()
2940
2941 @*/
PCFieldSplitGetDMSplits(PC pc,PetscBool * flg)2942 PetscErrorCode PCFieldSplitGetDMSplits(PC pc,PetscBool* flg)
2943 {
2944 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2945 PetscBool isfs;
2946 PetscErrorCode ierr;
2947
2948 PetscFunctionBegin;
2949 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
2950 PetscValidBoolPointer(flg,2);
2951 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr);
2952 if (isfs) {
2953 if (flg) *flg = jac->dm_splits;
2954 }
2955 PetscFunctionReturn(0);
2956 }
2957
2958 /*@
2959 PCFieldSplitGetDetectSaddlePoint - Returns flag indicating whether PCFieldSplit will attempt to automatically determine fields based on zero diagonal entries.
2960
2961 Logically Collective
2962
2963 Input Parameter:
2964 . pc - the preconditioner context
2965
2966 Output Parameter:
2967 . flg - boolean indicating whether to detect fields or not
2968
2969 Level: Intermediate
2970
2971 .seealso: PCFIELDSPLIT, PCFieldSplitSetDetectSaddlePoint()
2972
2973 @*/
PCFieldSplitGetDetectSaddlePoint(PC pc,PetscBool * flg)2974 PetscErrorCode PCFieldSplitGetDetectSaddlePoint(PC pc,PetscBool *flg)
2975 {
2976 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2977
2978 PetscFunctionBegin;
2979 *flg = jac->detect;
2980 PetscFunctionReturn(0);
2981 }
2982
2983 /*@
2984 PCFieldSplitSetDetectSaddlePoint - Sets flag indicating whether PCFieldSplit will attempt to automatically determine fields based on zero diagonal entries.
2985
2986 Logically Collective
2987
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()).
2990
2991 Input Parameter:
2992 . pc - the preconditioner context
2993
2994 Output Parameter:
2995 . flg - boolean indicating whether to detect fields or not
2996
2997 Options Database Key:
2998 . -pc_fieldsplit_detect_saddle_point
2999
3000 Level: Intermediate
3001
3002 .seealso: PCFIELDSPLIT, PCFieldSplitSetDetectSaddlePoint(), PCFieldSplitSetType(), PCFieldSplitSetSchurPre()
3003
3004 @*/
PCFieldSplitSetDetectSaddlePoint(PC pc,PetscBool flg)3005 PetscErrorCode PCFieldSplitSetDetectSaddlePoint(PC pc,PetscBool flg)
3006 {
3007 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
3008 PetscErrorCode ierr;
3009
3010 PetscFunctionBegin;
3011 jac->detect = flg;
3012 if (jac->detect) {
3013 ierr = PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);CHKERRQ(ierr);
3014 ierr = PCFieldSplitSetSchurPre(pc,PC_FIELDSPLIT_SCHUR_PRE_SELF,NULL);CHKERRQ(ierr);
3015 }
3016 PetscFunctionReturn(0);
3017 }
3018
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.
3023
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
3026
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
3029
3030 Level: intermediate
3031
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
3041
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_
3045
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.
3049
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().
3054
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.
3066
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.
3079
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.
3082
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.
3085
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".
3089
3090 References:
3091 A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier-Stokes equations
3092 Howard Elman, V.E. Howle, John Shadid, Robert Shuttleworth, Ray Tuminaro, Journal of Computational Physics 227 (2008) 1790--1808
3093 http://chess.cs.umd.edu/~elman/papers/tax.pdf
3094
3095 The Constrained Pressure Preconditioner (CPR) can be implemented using PCCOMPOSITE with PCGALERKIN. CPR first solves an R A P subsystem, updates the
3096 residual on all variables (PCCompositeSetType(pc,PC_COMPOSITE_MULTIPLICATIVE)), and then applies a simple ILU like preconditioner on all the variables.
3097
3098 The generalized Golub-Kahan bidiagonalization preconditioner (gkb) can be applied to symmetric 2x2 block matrices of the shape
3099 $ [ A00 A01]
3100 $ [ A01' 0 ]
3101 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'.
3102 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_.
3103
3104 [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013.
3105
3106 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCLSC,
3107 PCFieldSplitGetSubKSP(), PCFieldSplitSchurGetSubKSP(), PCFieldSplitSetFields(),
3108 PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre(), PCFieldSplitSetSchurFactType(),
3109 MatSchurComplementSetAinvType(), PCFieldSplitSetSchurScale(), PCFieldSplitSetDetectSaddlePoint()
3110 M*/
3111
PCCreate_FieldSplit(PC pc)3112 PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc)
3113 {
3114 PetscErrorCode ierr;
3115 PC_FieldSplit *jac;
3116
3117 PetscFunctionBegin;
3118 ierr = PetscNewLog(pc,&jac);CHKERRQ(ierr);
3119
3120 jac->bs = -1;
3121 jac->nsplits = 0;
3122 jac->type = PC_COMPOSITE_MULTIPLICATIVE;
3123 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */
3124 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL;
3125 jac->schurscale = -1.0;
3126 jac->dm_splits = PETSC_TRUE;
3127 jac->detect = PETSC_FALSE;
3128 jac->gkbtol = 1e-5;
3129 jac->gkbdelay = 5;
3130 jac->gkbnu = 1;
3131 jac->gkbmaxit = 100;
3132 jac->gkbmonitor = PETSC_FALSE;
3133
3134 pc->data = (void*)jac;
3135
3136 pc->ops->apply = PCApply_FieldSplit;
3137 pc->ops->applytranspose = PCApplyTranspose_FieldSplit;
3138 pc->ops->setup = PCSetUp_FieldSplit;
3139 pc->ops->reset = PCReset_FieldSplit;
3140 pc->ops->destroy = PCDestroy_FieldSplit;
3141 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit;
3142 pc->ops->view = PCView_FieldSplit;
3143 pc->ops->applyrichardson = NULL;
3144
3145 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurGetSubKSP_C",PCFieldSplitSchurGetSubKSP_FieldSplit);CHKERRQ(ierr);
3146 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr);
3147 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr);
3148 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr);
3149 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr);
3150 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr);
3151 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",PCFieldSplitRestrictIS_FieldSplit);CHKERRQ(ierr);
3152 PetscFunctionReturn(0);
3153 }
3154