Actual source code: fieldsplit.c
petsc-3.3-p7 2013-05-11
2: #include <petsc-private/pcimpl.h> /*I "petscpc.h" I*/
3: #include <petscdmcomposite.h> /*I "petscdmcomposite.h" I*/
5: const char *const PCFieldSplitSchurPreTypes[] = {"SELF","DIAG","USER","PCFieldSplitSchurPreType","PC_FIELDSPLIT_SCHUR_PRE_",0};
6: const char *const PCFieldSplitSchurFactTypes[] = {"DIAG","LOWER","UPPER","FULL","PCFieldSplitSchurFactType","PC_FIELDSPLIT_SCHUR_FACT_",0};
8: typedef struct _PC_FieldSplitLink *PC_FieldSplitLink;
9: struct _PC_FieldSplitLink {
10: KSP ksp;
11: Vec x,y;
12: char *splitname;
13: PetscInt nfields;
14: PetscInt *fields,*fields_col;
15: VecScatter sctx;
16: IS is,is_col;
17: PC_FieldSplitLink next,previous;
18: };
20: typedef struct {
21: PCCompositeType type;
22: PetscBool defaultsplit; /* Flag for a system with a set of 'k' scalar fields with the same layout (and bs = k) */
23: PetscBool splitdefined; /* Flag is set after the splits have been defined, to prevent more splits from being added */
24: PetscBool realdiagonal; /* Flag to use the diagonal blocks of mat preconditioned by pmat, instead of just pmat */
25: PetscInt bs; /* Block size for IS and Mat structures */
26: PetscInt nsplits; /* Number of field divisions defined */
27: Vec *x,*y,w1,w2;
28: Mat *mat; /* The diagonal block for each split */
29: Mat *pmat; /* The preconditioning diagonal block for each split */
30: Mat *Afield; /* The rows of the matrix associated with each split */
31: PetscBool issetup;
32: /* Only used when Schur complement preconditioning is used */
33: Mat B; /* The (0,1) block */
34: Mat C; /* The (1,0) block */
35: Mat schur; /* The Schur complement S = A11 - A10 A00^{-1} A01 */
36: Mat schur_user; /* User-provided preconditioning matrix for the Schur complement */
37: PCFieldSplitSchurPreType schurpre; /* Determines which preconditioning matrix is used for the Schur complement */
38: PCFieldSplitSchurFactType schurfactorization;
39: KSP kspschur; /* The solver for S */
40: PC_FieldSplitLink head;
41: PetscBool reset; /* indicates PCReset() has been last called on this object, hack */
42: PetscBool suboptionsset; /* Indicates that the KSPSetFromOptions() has been called on the sub-KSPs */
43: } PC_FieldSplit;
45: /*
46: Notes: there is no particular reason that pmat, x, and y are stored as arrays in PC_FieldSplit instead of
47: inside PC_FieldSplitLink, just historical. If you want to be able to add new fields after already using the
48: PC you could change this.
49: */
51: /* This helper is so that setting a user-provided preconditioning matrix is orthogonal to choosing to use it. This way the
52: * application-provided FormJacobian can provide this matrix without interfering with the user's (command-line) choices. */
53: static Mat FieldSplitSchurPre(PC_FieldSplit *jac)
54: {
55: switch (jac->schurpre) {
56: case PC_FIELDSPLIT_SCHUR_PRE_SELF: return jac->schur;
57: case PC_FIELDSPLIT_SCHUR_PRE_DIAG: return jac->pmat[1];
58: case PC_FIELDSPLIT_SCHUR_PRE_USER: /* Use a user-provided matrix if it is given, otherwise diagonal block */
59: default:
60: return jac->schur_user ? jac->schur_user : jac->pmat[1];
61: }
62: }
67: static PetscErrorCode PCView_FieldSplit(PC pc,PetscViewer viewer)
68: {
69: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
70: PetscErrorCode ierr;
71: PetscBool iascii;
72: PetscInt i,j;
73: PC_FieldSplitLink ilink = jac->head;
76: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
77: if (iascii) {
78: if (jac->bs > 0) {
79: PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs);
80: } else {
81: PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D\n",PCCompositeTypes[jac->type],jac->nsplits);
82: }
83: if (jac->realdiagonal) {
84: PetscViewerASCIIPrintf(viewer," using actual matrix for blocks rather than preconditioner matrix\n");
85: }
86: PetscViewerASCIIPrintf(viewer," Solver info for each split is in the following KSP objects:\n");
87: PetscViewerASCIIPushTab(viewer);
88: for (i=0; i<jac->nsplits; i++) {
89: if (ilink->fields) {
90: PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);
91: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
92: for (j=0; j<ilink->nfields; j++) {
93: if (j > 0) {
94: PetscViewerASCIIPrintf(viewer,",");
95: }
96: PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);
97: }
98: PetscViewerASCIIPrintf(viewer,"\n");
99: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
100: } else {
101: PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);
102: }
103: KSPView(ilink->ksp,viewer);
104: ilink = ilink->next;
105: }
106: PetscViewerASCIIPopTab(viewer);
107: } else {
108: SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for PCFieldSplit",((PetscObject)viewer)->type_name);
109: }
110: return(0);
111: }
115: static PetscErrorCode PCView_FieldSplit_Schur(PC pc,PetscViewer viewer)
116: {
117: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
118: PetscErrorCode ierr;
119: PetscBool iascii;
120: PetscInt i,j;
121: PC_FieldSplitLink ilink = jac->head;
122: KSP ksp;
125: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
126: if (iascii) {
127: if (jac->bs > 0) {
128: PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, blocksize = %D, factorization %s\n",jac->bs,PCFieldSplitSchurFactTypes[jac->schurfactorization]);
129: } else {
130: PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, factorization %s\n",PCFieldSplitSchurFactTypes[jac->schurfactorization]);
131: }
132: if (jac->realdiagonal) {
133: PetscViewerASCIIPrintf(viewer," using actual matrix for blocks rather than preconditioner matrix\n");
134: }
135: switch(jac->schurpre) {
136: case PC_FIELDSPLIT_SCHUR_PRE_SELF:
137: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from S itself\n");break;
138: case PC_FIELDSPLIT_SCHUR_PRE_DIAG:
139: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from the block diagonal part of A11\n");break;
140: case PC_FIELDSPLIT_SCHUR_PRE_USER:
141: if (jac->schur_user) {
142: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from user provided matrix\n");
143: } else {
144: PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from the block diagonal part of A11\n");
145: }
146: break;
147: default:
148: SETERRQ1(((PetscObject) pc)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid Schur preconditioning type: %d", jac->schurpre);
149: }
150: PetscViewerASCIIPrintf(viewer," Split info:\n");
151: PetscViewerASCIIPushTab(viewer);
152: for (i=0; i<jac->nsplits; i++) {
153: if (ilink->fields) {
154: PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);
155: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
156: for (j=0; j<ilink->nfields; j++) {
157: if (j > 0) {
158: PetscViewerASCIIPrintf(viewer,",");
159: }
160: PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);
161: }
162: PetscViewerASCIIPrintf(viewer,"\n");
163: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
164: } else {
165: PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);
166: }
167: ilink = ilink->next;
168: }
169: PetscViewerASCIIPrintf(viewer,"KSP solver for A00 block \n");
170: PetscViewerASCIIPushTab(viewer);
171: if (jac->schur) {
172: MatSchurComplementGetKSP(jac->schur,&ksp);
173: KSPView(ksp,viewer);
174: } else {
175: PetscViewerASCIIPrintf(viewer," not yet available\n");
176: }
177: PetscViewerASCIIPopTab(viewer);
178: PetscViewerASCIIPrintf(viewer,"KSP solver for S = A11 - A10 inv(A00) A01 \n");
179: PetscViewerASCIIPushTab(viewer);
180: if (jac->kspschur) {
181: KSPView(jac->kspschur,viewer);
182: } else {
183: PetscViewerASCIIPrintf(viewer," not yet available\n");
184: }
185: PetscViewerASCIIPopTab(viewer);
186: PetscViewerASCIIPopTab(viewer);
187: } else {
188: SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for PCFieldSplit",((PetscObject)viewer)->type_name);
189: }
190: return(0);
191: }
195: /* Precondition: jac->bs is set to a meaningful value */
196: static PetscErrorCode PCFieldSplitSetRuntimeSplits_Private(PC pc)
197: {
199: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
200: PetscInt i,nfields,*ifields,nfields_col,*ifields_col;
201: PetscBool flg,flg_col;
202: char optionname[128],splitname[8],optionname_col[128];
205: PetscMalloc(jac->bs*sizeof(PetscInt),&ifields);
206: PetscMalloc(jac->bs*sizeof(PetscInt),&ifields_col);
207: for (i=0,flg=PETSC_TRUE; ; i++) {
208: PetscSNPrintf(splitname,sizeof splitname,"%D",i);
209: PetscSNPrintf(optionname,sizeof optionname,"-pc_fieldsplit_%D_fields",i);
210: PetscSNPrintf(optionname_col,sizeof optionname_col,"-pc_fieldsplit_%D_fields_col",i);
211: nfields = jac->bs;
212: nfields_col = jac->bs;
213: PetscOptionsGetIntArray(((PetscObject)pc)->prefix,optionname,ifields,&nfields,&flg);
214: PetscOptionsGetIntArray(((PetscObject)pc)->prefix,optionname_col,ifields_col,&nfields_col,&flg_col);
215: if (!flg) break;
216: else if (flg && !flg_col) {
217: if (!nfields) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields");
218: PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields);
219: }
220: else {
221: if (!nfields || !nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields");
222: if (nfields != nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Number of row and column fields must match");
223: PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields_col);
224: }
225: }
226: if (i > 0) {
227: /* Makes command-line setting of splits take precedence over setting them in code.
228: Otherwise subsequent calls to PCFieldSplitSetIS() or PCFieldSplitSetFields() would
229: create new splits, which would probably not be what the user wanted. */
230: jac->splitdefined = PETSC_TRUE;
231: }
232: PetscFree(ifields);
233: PetscFree(ifields_col);
234: return(0);
235: }
239: static PetscErrorCode PCFieldSplitSetDefaults(PC pc)
240: {
241: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
242: PetscErrorCode ierr;
243: PC_FieldSplitLink ilink = jac->head;
244: PetscBool flg = PETSC_FALSE,stokes = PETSC_FALSE;
245: PetscInt i;
248: if (!ilink) {
249: PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,PETSC_NULL);
250: if (pc->dm && !stokes) {
251: PetscInt numFields, f;
252: char **fieldNames;
253: IS *fields;
254: DM *dms;
255: DMCreateFieldDecomposition(pc->dm, &numFields, &fieldNames, &fields, &dms);
256: for(f = 0; f < numFields; ++f) {
257: PCFieldSplitSetIS(pc, fieldNames[f], fields[f]);
258: PetscFree(fieldNames[f]);
259: ISDestroy(&fields[f]);
260: }
261: PetscFree(fieldNames);
262: PetscFree(fields);
263: if(dms) {
264: PetscInfo(pc,"Setting up physics based fieldsplit preconditioner using the embedded DM\n");
265: for (ilink=jac->head,i=0; ilink; ilink=ilink->next,i++) {
266: KSPSetDM(ilink->ksp, dms[i]);
267: KSPSetDMActive(ilink->ksp, PETSC_FALSE);
268: DMDestroy(&dms[i]);
269: }
270: PetscFree(dms);
271: }
272: } else {
273: if (jac->bs <= 0) {
274: if (pc->pmat) {
275: MatGetBlockSize(pc->pmat,&jac->bs);
276: } else {
277: jac->bs = 1;
278: }
279: }
281: PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&flg,PETSC_NULL);
282: if (stokes) {
283: IS zerodiags,rest;
284: PetscInt nmin,nmax;
286: MatGetOwnershipRange(pc->mat,&nmin,&nmax);
287: MatFindZeroDiagonals(pc->mat,&zerodiags);
288: ISComplement(zerodiags,nmin,nmax,&rest);
289: PCFieldSplitSetIS(pc,"0",rest);
290: PCFieldSplitSetIS(pc,"1",zerodiags);
291: ISDestroy(&zerodiags);
292: ISDestroy(&rest);
293: } else {
294: if (!flg) {
295: /* Allow user to set fields from command line, if bs was known at the time of PCSetFromOptions_FieldSplit()
296: then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */
297: PCFieldSplitSetRuntimeSplits_Private(pc);
298: if (jac->splitdefined) {PetscInfo(pc,"Splits defined using the options database\n");}
299: }
300: if (flg || !jac->splitdefined) {
301: PetscInfo(pc,"Using default splitting of fields\n");
302: for (i=0; i<jac->bs; i++) {
303: char splitname[8];
304: PetscSNPrintf(splitname,sizeof splitname,"%D",i);
305: PCFieldSplitSetFields(pc,splitname,1,&i,&i);
306: }
307: jac->defaultsplit = PETSC_TRUE;
308: }
309: }
310: }
311: } else if (jac->nsplits == 1) {
312: if (ilink->is) {
313: IS is2;
314: PetscInt nmin,nmax;
316: MatGetOwnershipRange(pc->mat,&nmin,&nmax);
317: ISComplement(ilink->is,nmin,nmax,&is2);
318: PCFieldSplitSetIS(pc,"1",is2);
319: ISDestroy(&is2);
320: } else SETERRQ(((PetscObject) pc)->comm,PETSC_ERR_SUP,"Must provide at least two sets of fields to PCFieldSplit()");
321: } else if (jac->reset) {
322: /* PCReset() has been called on this PC, ilink exists but all IS and Vec data structures in it must be rebuilt
323: This is basically the !ilink portion of code above copied from above and the allocation of the ilinks removed
324: since they already exist. This should be totally rewritten */
325: PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,PETSC_NULL);
326: if (pc->dm && !stokes) {
327: PetscBool dmcomposite;
328: PetscObjectTypeCompare((PetscObject)pc->dm,DMCOMPOSITE,&dmcomposite);
329: if (dmcomposite) {
330: PetscInt nDM;
331: IS *fields;
332: DM *dms;
333: PetscInfo(pc,"Setting up physics based fieldsplit preconditioner using the embedded DM\n");
334: DMCompositeGetNumberDM(pc->dm,&nDM);
335: DMCompositeGetGlobalISs(pc->dm,&fields);
336: PetscMalloc(nDM*sizeof(DM),&dms);
337: DMCompositeGetEntriesArray(pc->dm,dms);
338: for (i=0; i<nDM; i++) {
339: KSPSetDM(ilink->ksp,dms[i]);
340: KSPSetDMActive(ilink->ksp,PETSC_FALSE);
341: ilink->is = fields[i];
342: ilink = ilink->next;
343: }
344: PetscFree(fields);
345: PetscFree(dms);
346: }
347: } else {
348: PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&flg,PETSC_NULL);
349: if (stokes) {
350: IS zerodiags,rest;
351: PetscInt nmin,nmax;
353: MatGetOwnershipRange(pc->mat,&nmin,&nmax);
354: MatFindZeroDiagonals(pc->mat,&zerodiags);
355: ISComplement(zerodiags,nmin,nmax,&rest);
356: ISDestroy(&ilink->is);
357: ISDestroy(&ilink->next->is);
358: ilink->is = rest;
359: ilink->next->is = zerodiags;
360: } else SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Cases not yet handled when PCReset() was used");
361: }
362: }
364: if (jac->nsplits < 2) SETERRQ1(((PetscObject) pc)->comm,PETSC_ERR_PLIB,"Unhandled case, must have at least two fields, not %d", jac->nsplits);
365: return(0);
366: }
370: static PetscErrorCode PCSetUp_FieldSplit(PC pc)
371: {
372: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
373: PetscErrorCode ierr;
374: PC_FieldSplitLink ilink;
375: PetscInt i,nsplit;
376: MatStructure flag = pc->flag;
377: PetscBool sorted, sorted_col;
380: PCFieldSplitSetDefaults(pc);
381: nsplit = jac->nsplits;
382: ilink = jac->head;
384: /* get the matrices for each split */
385: if (!jac->issetup) {
386: PetscInt rstart,rend,nslots,bs;
388: jac->issetup = PETSC_TRUE;
390: /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */
391: if (jac->defaultsplit || !ilink->is) {
392: if (jac->bs <= 0) jac->bs = nsplit;
393: }
394: bs = jac->bs;
395: MatGetOwnershipRange(pc->pmat,&rstart,&rend);
396: nslots = (rend - rstart)/bs;
397: for (i=0; i<nsplit; i++) {
398: if (jac->defaultsplit) {
399: ISCreateStride(((PetscObject)pc)->comm,nslots,rstart+i,nsplit,&ilink->is);
400: ISDuplicate(ilink->is,&ilink->is_col);
401: } else if (!ilink->is) {
402: if (ilink->nfields > 1) {
403: PetscInt *ii,*jj,j,k,nfields = ilink->nfields,*fields = ilink->fields,*fields_col = ilink->fields_col;
404: PetscMalloc(ilink->nfields*nslots*sizeof(PetscInt),&ii);
405: PetscMalloc(ilink->nfields*nslots*sizeof(PetscInt),&jj);
406: for (j=0; j<nslots; j++) {
407: for (k=0; k<nfields; k++) {
408: ii[nfields*j + k] = rstart + bs*j + fields[k];
409: jj[nfields*j + k] = rstart + bs*j + fields_col[k];
410: }
411: }
412: ISCreateGeneral(((PetscObject)pc)->comm,nslots*nfields,ii,PETSC_OWN_POINTER,&ilink->is);
413: ISCreateGeneral(((PetscObject)pc)->comm,nslots*nfields,jj,PETSC_OWN_POINTER,&ilink->is_col);
414: } else {
415: ISCreateStride(((PetscObject)pc)->comm,nslots,rstart+ilink->fields[0],bs,&ilink->is);
416: ISCreateStride(((PetscObject)pc)->comm,nslots,rstart+ilink->fields_col[0],bs,&ilink->is_col);
417: }
418: }
419: ISSorted(ilink->is,&sorted);
420: if (ilink->is_col) { ISSorted(ilink->is_col,&sorted_col); }
421: if (!sorted || !sorted_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Fields must be sorted when creating split");
422: ilink = ilink->next;
423: }
424: }
426: ilink = jac->head;
427: if (!jac->pmat) {
428: Vec xtmp;
430: MatGetVecs(pc->pmat,&xtmp,PETSC_NULL);
431: PetscMalloc(nsplit*sizeof(Mat),&jac->pmat);
432: PetscMalloc2(nsplit,Vec,&jac->x,nsplit,Vec,&jac->y);
433: for (i=0; i<nsplit; i++) {
434: MatNullSpace sp;
436: /* Check for preconditioning matrix attached to IS */
437: PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject *) &jac->pmat[i]);
438: if (jac->pmat[i]) {
439: PetscObjectReference((PetscObject) jac->pmat[i]);
440: if (jac->type == PC_COMPOSITE_SCHUR) {
441: jac->schur_user = jac->pmat[i];
442: PetscObjectReference((PetscObject) jac->schur_user);
443: }
444: } else {
445: MatGetSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->pmat[i]);
446: }
447: /* create work vectors for each split */
448: MatGetVecs(jac->pmat[i],&jac->x[i],&jac->y[i]);
449: ilink->x = jac->x[i]; ilink->y = jac->y[i];
450: /* compute scatter contexts needed by multiplicative versions and non-default splits */
451: VecScatterCreate(xtmp,ilink->is,jac->x[i],PETSC_NULL,&ilink->sctx);
452: /* HACK: Check for the constant null space */
453: MatGetNullSpace(pc->pmat, &sp);
454: if (sp) {
455: MatNullSpace subsp;
456: Vec ftmp, gtmp;
457: PetscReal norm;
458: PetscInt N;
460: MatGetVecs(pc->pmat, >mp, PETSC_NULL);
461: MatGetVecs(jac->pmat[i], &ftmp, PETSC_NULL);
462: VecGetSize(ftmp, &N);
463: VecSet(ftmp, 1.0/N);
464: VecScatterBegin(ilink->sctx, ftmp, gtmp, INSERT_VALUES, SCATTER_REVERSE);
465: VecScatterEnd(ilink->sctx, ftmp, gtmp, INSERT_VALUES, SCATTER_REVERSE);
466: MatNullSpaceRemove(sp, gtmp, PETSC_NULL);
467: VecNorm(gtmp, NORM_2, &norm);
468: if (norm < 1.0e-10) {
469: MatNullSpaceCreate(((PetscObject)pc)->comm, PETSC_TRUE, 0, PETSC_NULL, &subsp);
470: MatSetNullSpace(jac->pmat[i], subsp);
471: MatNullSpaceDestroy(&subsp);
472: }
473: VecDestroy(&ftmp);
474: VecDestroy(>mp);
475: }
476: /* Check for null space attached to IS */
477: PetscObjectQuery((PetscObject) ilink->is, "nearnullspace", (PetscObject *) &sp);
478: if (sp) {
479: MatSetNearNullSpace(jac->pmat[i], sp);
480: }
481: ilink = ilink->next;
482: }
483: VecDestroy(&xtmp);
484: } else {
485: for (i=0; i<nsplit; i++) {
486: Mat pmat;
488: /* Check for preconditioning matrix attached to IS */
489: PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject *) &pmat);
490: if (!pmat) {
491: MatGetSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_REUSE_MATRIX,&jac->pmat[i]);
492: }
493: ilink = ilink->next;
494: }
495: }
496: if (jac->realdiagonal) {
497: ilink = jac->head;
498: if (!jac->mat) {
499: PetscMalloc(nsplit*sizeof(Mat),&jac->mat);
500: for (i=0; i<nsplit; i++) {
501: MatGetSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->mat[i]);
502: ilink = ilink->next;
503: }
504: } else {
505: for (i=0; i<nsplit; i++) {
506: if (jac->mat[i]) {MatGetSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_REUSE_MATRIX,&jac->mat[i]);}
507: ilink = ilink->next;
508: }
509: }
510: } else {
511: jac->mat = jac->pmat;
512: }
514: if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR) {
515: /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */
516: ilink = jac->head;
517: if (!jac->Afield) {
518: PetscMalloc(nsplit*sizeof(Mat),&jac->Afield);
519: for (i=0; i<nsplit; i++) {
520: MatGetSubMatrix(pc->mat,ilink->is,PETSC_NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);
521: ilink = ilink->next;
522: }
523: } else {
524: for (i=0; i<nsplit; i++) {
525: MatGetSubMatrix(pc->mat,ilink->is,PETSC_NULL,MAT_REUSE_MATRIX,&jac->Afield[i]);
526: ilink = ilink->next;
527: }
528: }
529: }
531: if (jac->type == PC_COMPOSITE_SCHUR) {
532: IS ccis;
533: PetscInt rstart,rend;
534: char lscname[256];
535: PetscObject LSC_L;
536: if (nsplit != 2) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_INCOMP,"To use Schur complement preconditioner you must have exactly 2 fields");
538: /* When extracting off-diagonal submatrices, we take complements from this range */
539: MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);
541: /* need to handle case when one is resetting up the preconditioner */
542: if (jac->schur) {
543: ilink = jac->head;
544: ISComplement(ilink->is_col,rstart,rend,&ccis);
545: MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);
546: ISDestroy(&ccis);
547: ilink = ilink->next;
548: ISComplement(ilink->is_col,rstart,rend,&ccis);
549: MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);
550: ISDestroy(&ccis);
551: MatSchurComplementUpdate(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1],pc->flag);
552: KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),pc->flag);
553: } else {
554: KSP ksp;
555: char schurprefix[256];
556: MatNullSpace sp;
558: /* extract the A01 and A10 matrices */
559: ilink = jac->head;
560: ISComplement(ilink->is_col,rstart,rend,&ccis);
561: MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
562: ISDestroy(&ccis);
563: ilink = ilink->next;
564: ISComplement(ilink->is_col,rstart,rend,&ccis);
565: MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
566: ISDestroy(&ccis);
567: /* Use mat[0] (diagonal block of the real matrix) preconditioned by pmat[0] */
568: MatCreateSchurComplement(jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1],&jac->schur);
569: MatGetNullSpace(jac->pmat[1], &sp);
570: /* FIXME: Attaching the A11-block's nullspace to S in lieu of a composable way to provide a null space for S */
571: if (sp) {MatSetNullSpace(jac->schur, sp);}
572: /* set tabbing, options prefix and DM of KSP inside the MatSchurComplement */
573: MatSchurComplementGetKSP(jac->schur,&ksp);
574: PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,2);
575: {
576: PC pcinner;
577: KSPGetPC(ksp,&pcinner);
578: PetscObjectIncrementTabLevel((PetscObject)pcinner,(PetscObject)pc,2);
579: }
580: PetscSNPrintf(schurprefix,sizeof schurprefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",jac->head->splitname);
581: KSPSetOptionsPrefix(ksp,schurprefix);
582: {
583: DM dminner;
584: KSPGetDM(jac->head->ksp,&dminner);
585: KSPSetDM(ksp,dminner);
586: KSPSetDMActive(ksp, PETSC_FALSE);
587: }
588: KSPSetFromOptions(ksp);
589: /* Need to call this everytime because new matrix is being created */
590: MatSetFromOptions(jac->schur);
591: MatSetUp(jac->schur);
592: /* Create and set up the KSP for the Schur complement; forward the second split's DM and set up tabbingn. */
593: KSPCreate(((PetscObject)pc)->comm,&jac->kspschur);
594: PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);
595: PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);
596: {
597: DM dmschur;
598: KSPGetDM(ilink->ksp,&dmschur);
599: KSPSetDM(jac->kspschur,dmschur);
600: KSPSetDMActive(jac->kspschur, PETSC_FALSE);
601: }
602: KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),DIFFERENT_NONZERO_PATTERN);
603: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) {
604: PC pcschur;
605: KSPGetPC(jac->kspschur,&pcschur);
606: PCSetType(pcschur,PCNONE);
607: /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */
608: }
609: PetscSNPrintf(schurprefix,sizeof schurprefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);
610: KSPSetOptionsPrefix(jac->kspschur,schurprefix);
611: /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
612: /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */
613: KSPSetFromOptions(jac->kspschur);
614: }
616: /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */
617: PetscSNPrintf(lscname,sizeof lscname,"%s_LSC_L",ilink->splitname);
618: PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);
619: if (!LSC_L) {PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);}
620: if (LSC_L) {PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L);}
621: PetscSNPrintf(lscname,sizeof lscname,"%s_LSC_Lp",ilink->splitname);
622: PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);
623: if (!LSC_L) {PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);}
624: if (LSC_L) {PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L);}
625: } else {
626: /* set up the individual splits' PCs */
627: i = 0;
628: ilink = jac->head;
629: while (ilink) {
630: KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i],flag);
631: /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
632: if (!jac->suboptionsset) {KSPSetFromOptions(ilink->ksp);}
633: i++;
634: ilink = ilink->next;
635: }
636: }
638: jac->suboptionsset = PETSC_TRUE;
639: return(0);
640: }
642: #define FieldSplitSplitSolveAdd(ilink,xx,yy) \
643: (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
644: VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
645: KSPSolve(ilink->ksp,ilink->x,ilink->y) || \
646: VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \
647: VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE))
651: static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y)
652: {
653: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
654: PetscErrorCode ierr;
655: KSP ksp;
656: PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
659: MatSchurComplementGetKSP(jac->schur,&ksp);
661: switch (jac->schurfactorization) {
662: case PC_FIELDSPLIT_SCHUR_FACT_DIAG:
663: /* [A00 0; 0 -S], positive definite, suitable for MINRES */
664: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
665: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
666: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
667: KSPSolve(ksp,ilinkA->x,ilinkA->y);
668: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
669: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
670: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
671: VecScale(ilinkD->y,-1.);
672: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
673: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
674: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
675: break;
676: case PC_FIELDSPLIT_SCHUR_FACT_LOWER:
677: /* [A00 0; A10 S], suitable for left preconditioning */
678: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
679: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
680: KSPSolve(ksp,ilinkA->x,ilinkA->y);
681: MatMult(jac->C,ilinkA->y,ilinkD->x);
682: VecScale(ilinkD->x,-1.);
683: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
684: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
685: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
686: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
687: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
688: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
689: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
690: break;
691: case PC_FIELDSPLIT_SCHUR_FACT_UPPER:
692: /* [A00 A01; 0 S], suitable for right preconditioning */
693: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
694: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
695: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
696: MatMult(jac->B,ilinkD->y,ilinkA->x);
697: VecScale(ilinkA->x,-1.);
698: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);
699: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
700: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);
701: KSPSolve(ksp,ilinkA->x,ilinkA->y);
702: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
703: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
704: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
705: break;
706: case PC_FIELDSPLIT_SCHUR_FACT_FULL:
707: /* [1 0; A10 A00^{-1} 1] [A00 0; 0 S] [1 A00^{-1}A01; 0 1], an exact solve if applied exactly, needs one extra solve with A */
708: VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
709: VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
710: KSPSolve(ksp,ilinkA->x,ilinkA->y);
711: MatMult(jac->C,ilinkA->y,ilinkD->x);
712: VecScale(ilinkD->x,-1.0);
713: VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
714: VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
716: KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
717: VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
718: VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
720: MatMult(jac->B,ilinkD->y,ilinkA->y);
721: VecAXPY(ilinkA->x,-1.0,ilinkA->y);
722: KSPSolve(ksp,ilinkA->x,ilinkA->y);
723: VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
724: VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
725: }
726: return(0);
727: }
731: static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y)
732: {
733: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
734: PetscErrorCode ierr;
735: PC_FieldSplitLink ilink = jac->head;
736: PetscInt cnt,bs;
739: CHKMEMQ;
740: if (jac->type == PC_COMPOSITE_ADDITIVE) {
741: if (jac->defaultsplit) {
742: VecGetBlockSize(x,&bs);
743: if (jac->bs > 0 && bs != jac->bs) SETERRQ2(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Blocksize of x vector %D does not match fieldsplit blocksize %D",bs,jac->bs);
744: VecGetBlockSize(y,&bs);
745: if (jac->bs > 0 && bs != jac->bs) SETERRQ2(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Blocksize of y vector %D does not match fieldsplit blocksize %D",bs,jac->bs);
746: VecStrideGatherAll(x,jac->x,INSERT_VALUES);
747: while (ilink) {
748: KSPSolve(ilink->ksp,ilink->x,ilink->y);
749: ilink = ilink->next;
750: }
751: VecStrideScatterAll(jac->y,y,INSERT_VALUES);
752: } else {
753: VecSet(y,0.0);
754: while (ilink) {
755: FieldSplitSplitSolveAdd(ilink,x,y);
756: ilink = ilink->next;
757: }
758: }
759: } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
760: if (!jac->w1) {
761: VecDuplicate(x,&jac->w1);
762: VecDuplicate(x,&jac->w2);
763: }
764: VecSet(y,0.0);
765: FieldSplitSplitSolveAdd(ilink,x,y);
766: cnt = 1;
767: while (ilink->next) {
768: ilink = ilink->next;
769: /* compute the residual only over the part of the vector needed */
770: MatMult(jac->Afield[cnt++],y,ilink->x);
771: VecScale(ilink->x,-1.0);
772: VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
773: VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
774: KSPSolve(ilink->ksp,ilink->x,ilink->y);
775: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
776: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
777: }
778: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
779: cnt -= 2;
780: while (ilink->previous) {
781: ilink = ilink->previous;
782: /* compute the residual only over the part of the vector needed */
783: MatMult(jac->Afield[cnt--],y,ilink->x);
784: VecScale(ilink->x,-1.0);
785: VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
786: VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
787: KSPSolve(ilink->ksp,ilink->x,ilink->y);
788: VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
789: VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
790: }
791: }
792: } else SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type);
793: CHKMEMQ;
794: return(0);
795: }
797: #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \
798: (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
799: VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
800: KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \
801: VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \
802: VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE))
806: static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y)
807: {
808: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
809: PetscErrorCode ierr;
810: PC_FieldSplitLink ilink = jac->head;
811: PetscInt bs;
814: CHKMEMQ;
815: if (jac->type == PC_COMPOSITE_ADDITIVE) {
816: if (jac->defaultsplit) {
817: VecGetBlockSize(x,&bs);
818: if (jac->bs > 0 && bs != jac->bs) SETERRQ2(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Blocksize of x vector %D does not match fieldsplit blocksize %D",bs,jac->bs);
819: VecGetBlockSize(y,&bs);
820: if (jac->bs > 0 && bs != jac->bs) SETERRQ2(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Blocksize of y vector %D does not match fieldsplit blocksize %D",bs,jac->bs);
821: VecStrideGatherAll(x,jac->x,INSERT_VALUES);
822: while (ilink) {
823: KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);
824: ilink = ilink->next;
825: }
826: VecStrideScatterAll(jac->y,y,INSERT_VALUES);
827: } else {
828: VecSet(y,0.0);
829: while (ilink) {
830: FieldSplitSplitSolveAddTranspose(ilink,x,y);
831: ilink = ilink->next;
832: }
833: }
834: } else {
835: if (!jac->w1) {
836: VecDuplicate(x,&jac->w1);
837: VecDuplicate(x,&jac->w2);
838: }
839: VecSet(y,0.0);
840: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
841: FieldSplitSplitSolveAddTranspose(ilink,x,y);
842: while (ilink->next) {
843: ilink = ilink->next;
844: MatMultTranspose(pc->mat,y,jac->w1);
845: VecWAXPY(jac->w2,-1.0,jac->w1,x);
846: FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
847: }
848: while (ilink->previous) {
849: ilink = ilink->previous;
850: MatMultTranspose(pc->mat,y,jac->w1);
851: VecWAXPY(jac->w2,-1.0,jac->w1,x);
852: FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
853: }
854: } else {
855: while (ilink->next) { /* get to last entry in linked list */
856: ilink = ilink->next;
857: }
858: FieldSplitSplitSolveAddTranspose(ilink,x,y);
859: while (ilink->previous) {
860: ilink = ilink->previous;
861: MatMultTranspose(pc->mat,y,jac->w1);
862: VecWAXPY(jac->w2,-1.0,jac->w1,x);
863: FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
864: }
865: }
866: }
867: CHKMEMQ;
868: return(0);
869: }
873: static PetscErrorCode PCReset_FieldSplit(PC pc)
874: {
875: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
876: PetscErrorCode ierr;
877: PC_FieldSplitLink ilink = jac->head,next;
880: while (ilink) {
881: KSPReset(ilink->ksp);
882: VecDestroy(&ilink->x);
883: VecDestroy(&ilink->y);
884: VecScatterDestroy(&ilink->sctx);
885: ISDestroy(&ilink->is);
886: ISDestroy(&ilink->is_col);
887: next = ilink->next;
888: ilink = next;
889: }
890: PetscFree2(jac->x,jac->y);
891: if (jac->mat && jac->mat != jac->pmat) {
892: MatDestroyMatrices(jac->nsplits,&jac->mat);
893: } else if (jac->mat) {
894: jac->mat = PETSC_NULL;
895: }
896: if (jac->pmat) {MatDestroyMatrices(jac->nsplits,&jac->pmat);}
897: if (jac->Afield) {MatDestroyMatrices(jac->nsplits,&jac->Afield);}
898: VecDestroy(&jac->w1);
899: VecDestroy(&jac->w2);
900: MatDestroy(&jac->schur);
901: MatDestroy(&jac->schur_user);
902: KSPDestroy(&jac->kspschur);
903: MatDestroy(&jac->B);
904: MatDestroy(&jac->C);
905: jac->reset = PETSC_TRUE;
906: return(0);
907: }
911: static PetscErrorCode PCDestroy_FieldSplit(PC pc)
912: {
913: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
914: PetscErrorCode ierr;
915: PC_FieldSplitLink ilink = jac->head,next;
918: PCReset_FieldSplit(pc);
919: while (ilink) {
920: KSPDestroy(&ilink->ksp);
921: next = ilink->next;
922: PetscFree(ilink->splitname);
923: PetscFree(ilink->fields);
924: PetscFree(ilink->fields_col);
925: PetscFree(ilink);
926: ilink = next;
927: }
928: PetscFree2(jac->x,jac->y);
929: PetscFree(pc->data);
930: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","",PETSC_NULL);
931: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","",PETSC_NULL);
932: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetIS_C","",PETSC_NULL);
933: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetType_C","",PETSC_NULL);
934: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetBlockSize_C","",PETSC_NULL);
935: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","",PETSC_NULL);
936: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetSchurFactType_C","",PETSC_NULL);
937: return(0);
938: }
942: static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc)
943: {
944: PetscErrorCode ierr;
945: PetscInt bs;
946: PetscBool flg,stokes = PETSC_FALSE;
947: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
948: PCCompositeType ctype;
949: DM ddm;
950: char ddm_name[1025];
953: PetscOptionsHead("FieldSplit options");
954: if(pc->dm) {
955: /* Allow the user to request a decomposition DM by name */
956: PetscStrncpy(ddm_name, "", 1024);
957: PetscOptionsString("-pc_fieldsplit_decomposition", "Name of the DM defining the composition", "PCSetDM", ddm_name, ddm_name,1024,&flg);
958: if(flg) {
959: DMCreateFieldDecompositionDM(pc->dm, ddm_name, &ddm);
960: if(!ddm) {
961: SETERRQ1(((PetscObject)pc)->comm, PETSC_ERR_ARG_WRONGSTATE, "Uknown field decomposition name %s", ddm_name);
962: }
963: PetscInfo(pc,"Using field decomposition DM defined using options database\n");
964: PCSetDM(pc,ddm);
965: }
966: }
967: PetscOptionsBool("-pc_fieldsplit_real_diagonal","Use diagonal blocks of the operator","PCFieldSplitSetRealDiagonal",jac->realdiagonal,&jac->realdiagonal,PETSC_NULL);
968: PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);
969: if (flg) {
970: PCFieldSplitSetBlockSize(pc,bs);
971: }
973: PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,PETSC_NULL);
974: if (stokes) {
975: PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);
976: jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF;
977: }
979: PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);
980: if (flg) {
981: PCFieldSplitSetType(pc,ctype);
982: }
983: /* Only setup fields once */
984: if ((jac->bs > 0) && (jac->nsplits == 0)) {
985: /* only allow user to set fields from command line if bs is already known.
986: otherwise user can set them in PCFieldSplitSetDefaults() */
987: PCFieldSplitSetRuntimeSplits_Private(pc);
988: if (jac->splitdefined) {PetscInfo(pc,"Splits defined using the options database\n");}
989: }
990: if (jac->type == PC_COMPOSITE_SCHUR) {
991: PetscOptionsGetEnum(((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);
992: if (flg) {PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");}
993: PetscOptionsEnum("-pc_fieldsplit_schur_fact_type","Which off-diagonal parts of the block factorization to use","PCFieldSplitSetSchurFactType",PCFieldSplitSchurFactTypes,(PetscEnum)jac->schurfactorization,(PetscEnum*)&jac->schurfactorization,PETSC_NULL);
994: PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSchurPrecondition",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,PETSC_NULL);
995: }
996: PetscOptionsTail();
997: return(0);
998: }
1000: /*------------------------------------------------------------------------------------*/
1002: EXTERN_C_BEGIN
1005: PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1006: {
1007: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1008: PetscErrorCode ierr;
1009: PC_FieldSplitLink ilink,next = jac->head;
1010: char prefix[128];
1011: PetscInt i;
1014: if (jac->splitdefined) {
1015: PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1016: return(0);
1017: }
1018: for (i=0; i<n; i++) {
1019: if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs);
1020: if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]);
1021: }
1022: PetscNew(struct _PC_FieldSplitLink,&ilink);
1023: if (splitname) {
1024: PetscStrallocpy(splitname,&ilink->splitname);
1025: } else {
1026: PetscMalloc(3*sizeof(char),&ilink->splitname);
1027: PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);
1028: }
1029: PetscMalloc(n*sizeof(PetscInt),&ilink->fields);
1030: PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));
1031: PetscMalloc(n*sizeof(PetscInt),&ilink->fields_col);
1032: PetscMemcpy(ilink->fields_col,fields_col,n*sizeof(PetscInt));
1033: ilink->nfields = n;
1034: ilink->next = PETSC_NULL;
1035: KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);
1036: PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1037: KSPSetType(ilink->ksp,KSPPREONLY);
1038: PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);
1040: PetscSNPrintf(prefix,sizeof prefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);
1041: KSPSetOptionsPrefix(ilink->ksp,prefix);
1043: if (!next) {
1044: jac->head = ilink;
1045: ilink->previous = PETSC_NULL;
1046: } else {
1047: while (next->next) {
1048: next = next->next;
1049: }
1050: next->next = ilink;
1051: ilink->previous = next;
1052: }
1053: jac->nsplits++;
1054: return(0);
1055: }
1056: EXTERN_C_END
1058: EXTERN_C_BEGIN
1061: PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp)
1062: {
1063: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1067: PetscMalloc(jac->nsplits*sizeof(KSP),subksp);
1068: MatSchurComplementGetKSP(jac->schur,*subksp);
1069: (*subksp)[1] = jac->kspschur;
1070: if (n) *n = jac->nsplits;
1071: return(0);
1072: }
1073: EXTERN_C_END
1075: EXTERN_C_BEGIN
1078: PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp)
1079: {
1080: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1081: PetscErrorCode ierr;
1082: PetscInt cnt = 0;
1083: PC_FieldSplitLink ilink = jac->head;
1086: PetscMalloc(jac->nsplits*sizeof(KSP),subksp);
1087: while (ilink) {
1088: (*subksp)[cnt++] = ilink->ksp;
1089: ilink = ilink->next;
1090: }
1091: 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);
1092: if (n) *n = jac->nsplits;
1093: return(0);
1094: }
1095: EXTERN_C_END
1097: EXTERN_C_BEGIN
1100: PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is)
1101: {
1102: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1103: PetscErrorCode ierr;
1104: PC_FieldSplitLink ilink, next = jac->head;
1105: char prefix[128];
1108: if (jac->splitdefined) {
1109: PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1110: return(0);
1111: }
1112: PetscNew(struct _PC_FieldSplitLink,&ilink);
1113: if (splitname) {
1114: PetscStrallocpy(splitname,&ilink->splitname);
1115: } else {
1116: PetscMalloc(8*sizeof(char),&ilink->splitname);
1117: PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);
1118: }
1119: PetscObjectReference((PetscObject)is);
1120: ISDestroy(&ilink->is);
1121: ilink->is = is;
1122: PetscObjectReference((PetscObject)is);
1123: ISDestroy(&ilink->is_col);
1124: ilink->is_col = is;
1125: ilink->next = PETSC_NULL;
1126: KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);
1127: PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1128: KSPSetType(ilink->ksp,KSPPREONLY);
1129: PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);
1131: PetscSNPrintf(prefix,sizeof prefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);
1132: KSPSetOptionsPrefix(ilink->ksp,prefix);
1134: if (!next) {
1135: jac->head = ilink;
1136: ilink->previous = PETSC_NULL;
1137: } else {
1138: while (next->next) {
1139: next = next->next;
1140: }
1141: next->next = ilink;
1142: ilink->previous = next;
1143: }
1144: jac->nsplits++;
1146: return(0);
1147: }
1148: EXTERN_C_END
1152: /*@
1153: PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner
1155: Logically Collective on PC
1157: Input Parameters:
1158: + pc - the preconditioner context
1159: . splitname - name of this split, if PETSC_NULL the number of the split is used
1160: . n - the number of fields in this split
1161: - fields - the fields in this split
1163: Level: intermediate
1165: Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field.
1167: The PCFieldSplitSetFields() is for defining fields as a strided blocks. For example, if the block
1168: 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
1169: 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x....
1170: where the numbered entries indicate what is in the field.
1172: This function is called once per split (it creates a new split each time). Solve options
1173: for this split will be available under the prefix -fieldsplit_SPLITNAME_.
1175: Developer Note: This routine does not actually create the IS representing the split, that is delayed
1176: until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be
1177: available when this routine is called.
1179: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS()
1181: @*/
1182: PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1183: {
1189: if (n < 1) SETERRQ2(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname);
1191: PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt *,const PetscInt *),(pc,splitname,n,fields,fields_col));
1192: return(0);
1193: }
1197: /*@
1198: PCFieldSplitSetIS - Sets the exact elements for field
1200: Logically Collective on PC
1202: Input Parameters:
1203: + pc - the preconditioner context
1204: . splitname - name of this split, if PETSC_NULL the number of the split is used
1205: - is - the index set that defines the vector elements in this field
1208: Notes:
1209: Use PCFieldSplitSetFields(), for fields defined by strided types.
1211: This function is called once per split (it creates a new split each time). Solve options
1212: for this split will be available under the prefix -fieldsplit_SPLITNAME_.
1214: Level: intermediate
1216: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize()
1218: @*/
1219: PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is)
1220: {
1227: PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));
1228: return(0);
1229: }
1233: /*@
1234: PCFieldSplitGetIS - Retrieves the elements for a field as an IS
1236: Logically Collective on PC
1238: Input Parameters:
1239: + pc - the preconditioner context
1240: - splitname - name of this split
1242: Output Parameter:
1243: - is - the index set that defines the vector elements in this field, or PETSC_NULL if the field is not found
1245: Level: intermediate
1247: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS()
1249: @*/
1250: PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is)
1251: {
1258: {
1259: PC_FieldSplit *jac = (PC_FieldSplit *) pc->data;
1260: PC_FieldSplitLink ilink = jac->head;
1261: PetscBool found;
1263: *is = PETSC_NULL;
1264: while(ilink) {
1265: PetscStrcmp(ilink->splitname, splitname, &found);
1266: if (found) {
1267: *is = ilink->is;
1268: break;
1269: }
1270: ilink = ilink->next;
1271: }
1272: }
1273: return(0);
1274: }
1278: /*@
1279: PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the
1280: fieldsplit preconditioner. If not set the matrix block size is used.
1282: Logically Collective on PC
1284: Input Parameters:
1285: + pc - the preconditioner context
1286: - bs - the block size
1288: Level: intermediate
1290: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields()
1292: @*/
1293: PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs)
1294: {
1300: PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));
1301: return(0);
1302: }
1306: /*@C
1307: PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits
1308:
1309: Collective on KSP
1311: Input Parameter:
1312: . pc - the preconditioner context
1314: Output Parameters:
1315: + n - the number of splits
1316: - pc - the array of KSP contexts
1318: Note:
1319: After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user
1320: (not the KSP just the array that contains them).
1322: You must call KSPSetUp() before calling PCFieldSplitGetSubKSP().
1324: Level: advanced
1326: .seealso: PCFIELDSPLIT
1327: @*/
1328: PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
1329: {
1335: PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));
1336: return(0);
1337: }
1341: /*@
1342: PCFieldSplitSchurPrecondition - Indicates if the Schur complement is preconditioned by a preconditioner constructed by the
1343: A11 matrix. Otherwise no preconditioner is used.
1345: Collective on PC
1347: Input Parameters:
1348: + pc - the preconditioner context
1349: . ptype - which matrix to use for preconditioning the Schur complement, PC_FIELDSPLIT_SCHUR_PRE_DIAG (diag) is default
1350: - userpre - matrix to use for preconditioning, or PETSC_NULL
1352: Options Database:
1353: . -pc_fieldsplit_schur_precondition <self,user,diag> default is diag
1355: Notes:
1356: $ If ptype is
1357: $ user then the preconditioner for the Schur complement is generated by the provided matrix (pre argument
1358: $ to this function).
1359: $ diag then the preconditioner for the Schur complement is generated by the block diagonal part of the original
1360: $ matrix associated with the Schur complement (i.e. A11)
1361: $ self the preconditioner for the Schur complement is generated from the Schur complement matrix itself:
1362: $ The only preconditioner that currently works directly with the Schur complement matrix object is the PCLSC
1363: $ preconditioner
1365: When solving a saddle point problem, where the A11 block is identically zero, using diag as the ptype only makes sense
1366: with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and
1367: -fieldsplit_1_pc_type lsc which uses the least squares commutator compute a preconditioner for the Schur complement.
1368:
1369: Developer Notes: This is a terrible name, gives no good indication of what the function does and should also have Set in
1370: the name since it sets a proceedure to use.
1371:
1372: Level: intermediate
1374: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC
1376: @*/
1377: PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
1378: {
1383: PetscTryMethod(pc,"PCFieldSplitSchurPrecondition_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));
1384: return(0);
1385: }
1387: EXTERN_C_BEGIN
1390: PetscErrorCode PCFieldSplitSchurPrecondition_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
1391: {
1392: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1393: PetscErrorCode ierr;
1396: jac->schurpre = ptype;
1397: if (pre) {
1398: MatDestroy(&jac->schur_user);
1399: jac->schur_user = pre;
1400: PetscObjectReference((PetscObject)jac->schur_user);
1401: }
1402: return(0);
1403: }
1404: EXTERN_C_END
1408: /*@
1409: PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain
1411: Collective on PC
1413: Input Parameters:
1414: + pc - the preconditioner context
1415: - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default
1417: Options Database:
1418: . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full
1421: Level: intermediate
1423: Notes:
1424: The FULL factorization is
1426: $ (A B) = (1 0) (A 0) (1 Ainv*B)
1427: $ (C D) (C*Ainv 1) (0 S) (0 1 )
1429: where S = D - 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,
1430: 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).
1432: If applied exactly, FULL factorization is a direct solver. The preconditioned operator with LOWER or UPPER has all eigenvalues equal to 1 and minimal polynomial
1433: of degree 2, so KSPGMRES converges in 2 iterations. If the iteration count is very low, consider using KSPFGMRES or KSPGCR which can use one less preconditioner
1434: application in this case. Note that the preconditioned operator may be highly non-normal, so such fast convergence may not be observed in practice. With DIAG,
1435: the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so KSPGMRES converges in at most 4 iterations.
1437: For symmetric problems in which A is positive definite and S is negative definite, DIAG can be used with KSPMINRES. Note that a flexible method like KSPFGMRES
1438: or KSPGCR must be used if the fieldsplit preconditioner is nonlinear (e.g. a few iterations of a Krylov method is used inside a split).
1440: References:
1441: Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000) pp. 1969-1972.
1443: Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001), pp. 1050-1051.
1445: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType
1446: @*/
1447: PetscErrorCode PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype)
1448: {
1453: PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));
1454: return(0);
1455: }
1460: {
1461: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1464: jac->schurfactorization = ftype;
1465: return(0);
1466: }
1470: /*@C
1471: PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement
1472:
1473: Collective on KSP
1475: Input Parameter:
1476: . pc - the preconditioner context
1478: Output Parameters:
1479: + A00 - the (0,0) block
1480: . A01 - the (0,1) block
1481: . A10 - the (1,0) block
1482: - A11 - the (1,1) block
1484: Level: advanced
1486: .seealso: PCFIELDSPLIT
1487: @*/
1488: PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11)
1489: {
1490: PC_FieldSplit *jac = (PC_FieldSplit *) pc->data;
1494: if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach.");
1495: if (A00) *A00 = jac->pmat[0];
1496: if (A01) *A01 = jac->B;
1497: if (A10) *A10 = jac->C;
1498: if (A11) *A11 = jac->pmat[1];
1499: return(0);
1500: }
1502: EXTERN_C_BEGIN
1505: PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type)
1506: {
1507: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1511: jac->type = type;
1512: if (type == PC_COMPOSITE_SCHUR) {
1513: pc->ops->apply = PCApply_FieldSplit_Schur;
1514: pc->ops->view = PCView_FieldSplit_Schur;
1515: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit_Schur",PCFieldSplitGetSubKSP_FieldSplit_Schur);
1516: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","PCFieldSplitSchurPrecondition_FieldSplit",PCFieldSplitSchurPrecondition_FieldSplit);
1517: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetSchurFactType_C","PCFieldSplitSetSchurFactType_FieldSplit",PCFieldSplitSetSchurFactType_FieldSplit);
1519: } else {
1520: pc->ops->apply = PCApply_FieldSplit;
1521: pc->ops->view = PCView_FieldSplit;
1522: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit",PCFieldSplitGetSubKSP_FieldSplit);
1523: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","",0);
1524: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetSchurFactType_C","",0);
1525: }
1526: return(0);
1527: }
1528: EXTERN_C_END
1530: EXTERN_C_BEGIN
1533: PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs)
1534: {
1535: PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1538: if (bs < 1) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs);
1539: if (jac->bs > 0 && jac->bs != bs) SETERRQ2(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Cannot change fieldsplit blocksize from %D to %D after it has been set",jac->bs,bs);
1540: jac->bs = bs;
1541: return(0);
1542: }
1543: EXTERN_C_END
1547: /*@
1548: PCFieldSplitSetType - Sets the type of fieldsplit preconditioner.
1549:
1550: Collective on PC
1552: Input Parameter:
1553: . pc - the preconditioner context
1554: . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR
1556: Options Database Key:
1557: . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type
1559: Level: Intermediate
1561: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
1563: .seealso: PCCompositeSetType()
1565: @*/
1566: PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type)
1567: {
1572: PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));
1573: return(0);
1574: }
1578: /*@
1579: PCFieldSplitGetType - Gets the type of fieldsplit preconditioner.
1581: Not collective
1583: Input Parameter:
1584: . pc - the preconditioner context
1586: Output Parameter:
1587: . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR
1589: Level: Intermediate
1591: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
1592: .seealso: PCCompositeSetType()
1593: @*/
1594: PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type)
1595: {
1596: PC_FieldSplit *jac = (PC_FieldSplit *) pc->data;
1601: *type = jac->type;
1602: return(0);
1603: }
1605: /* -------------------------------------------------------------------------------------*/
1606: /*MC
1607: PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual
1608: fields or groups of fields. See the users manual section "Solving Block Matrices" for more details.
1610: To set options on the solvers for each block append -fieldsplit_ to all the PC
1611: options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1
1612:
1613: To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP()
1614: and set the options directly on the resulting KSP object
1616: Level: intermediate
1618: Options Database Keys:
1619: + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split
1620: . -pc_fieldsplit_default - automatically add any fields to additional splits that have not
1621: been supplied explicitly by -pc_fieldsplit_%d_fields
1622: . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields)
1623: . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur> - type of relaxation or factorization splitting
1624: . -pc_fieldsplit_schur_precondition <self,user,diag> - default is diag
1625: . -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero or negative diagonal and uses Schur complement with no preconditioner as the solver
1627: - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_
1628: for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields
1630: Notes:
1631: Use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS()
1632: to define a field by an arbitrary collection of entries.
1634: If no fields are set the default is used. The fields are defined by entries strided by bs,
1635: beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(),
1636: if this is not called the block size defaults to the blocksize of the second matrix passed
1637: to KSPSetOperators()/PCSetOperators().
1639: $ For the Schur complement preconditioner if J = ( A00 A01 )
1640: $ ( A10 A11 )
1641: $ the preconditioner using full factorization is
1642: $ ( I -A10 ksp(A00) ) ( inv(A00) 0 ) ( I 0 )
1643: $ ( 0 I ) ( 0 ksp(S) ) ( -A10 ksp(A00) I )
1644: where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. The action of
1645: ksp(S) is computed using the KSP solver with prefix -fieldsplit_splitname_ (where splitname was given
1646: in providing the SECOND split or 1 if not give). For PCFieldSplitGetKSP() when field number is 0,
1647: it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default
1648: A11 is used to construct a preconditioner for S, use PCFieldSplitSchurPrecondition() to turn on or off this
1649: option. You can use the preconditioner PCLSC to precondition the Schur complement with -fieldsplit_1_pc_type lsc. The
1650: factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above,
1651: diag gives
1652: $ ( inv(A00) 0 )
1653: $ ( 0 -ksp(S) )
1654: note that slightly counter intuitively there is a negative in front of the ksp(S) so that the preconditioner is positive definite. The lower factorization is the inverse of
1655: $ ( A00 0 )
1656: $ ( A10 S )
1657: where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of
1658: $ ( A00 A01 )
1659: $ ( 0 S )
1660: where again the inverses of A00 and S are applied using KSPs.
1662: If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS
1663: is used automatically for a second block.
1665: The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1.
1666: Generally it should be used with the AIJ format.
1668: The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see,
1669: for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT
1670: inside a smoother resulting in "Distributive Smoothers".
1672: Concepts: physics based preconditioners, block preconditioners
1674: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC,
1675: PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSchurPrecondition()
1676: M*/
1678: EXTERN_C_BEGIN
1681: PetscErrorCode PCCreate_FieldSplit(PC pc)
1682: {
1684: PC_FieldSplit *jac;
1687: PetscNewLog(pc,PC_FieldSplit,&jac);
1688: jac->bs = -1;
1689: jac->nsplits = 0;
1690: jac->type = PC_COMPOSITE_MULTIPLICATIVE;
1691: jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */
1692: jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL;
1694: pc->data = (void*)jac;
1696: pc->ops->apply = PCApply_FieldSplit;
1697: pc->ops->applytranspose = PCApplyTranspose_FieldSplit;
1698: pc->ops->setup = PCSetUp_FieldSplit;
1699: pc->ops->reset = PCReset_FieldSplit;
1700: pc->ops->destroy = PCDestroy_FieldSplit;
1701: pc->ops->setfromoptions = PCSetFromOptions_FieldSplit;
1702: pc->ops->view = PCView_FieldSplit;
1703: pc->ops->applyrichardson = 0;
1705: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit",
1706: PCFieldSplitGetSubKSP_FieldSplit);
1707: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","PCFieldSplitSetFields_FieldSplit",
1708: PCFieldSplitSetFields_FieldSplit);
1709: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetIS_C","PCFieldSplitSetIS_FieldSplit",
1710: PCFieldSplitSetIS_FieldSplit);
1711: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetType_C","PCFieldSplitSetType_FieldSplit",
1712: PCFieldSplitSetType_FieldSplit);
1713: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetBlockSize_C","PCFieldSplitSetBlockSize_FieldSplit",
1714: PCFieldSplitSetBlockSize_FieldSplit);
1715: return(0);
1716: }
1717: EXTERN_C_END