Actual source code: bfbt.c
1: #include "private/pcimpl.h" /*I "petscpc.h" I*/
3: typedef struct {
4: Mat K; /* K, the (0,0) block, is [M x M] */
5: Mat G; /* G, the (0,1) block, is [M x N] */
6: Mat M; /* The mass matrix corresponding to K */
7: Vec inv_diag_M;
8: Mat GtG;
9: Vec s, t, X; /* s is [M], t is [N], X is [M] */
10: KSP ksp;
11: PetscTruth form_GtG; /* Don't allow anything else yet */
12: PetscTruth scaled; /* Use K for M */
13: } PC_BFBt;
18: /* K and G must not be PETSC_NULL, M can be PETSC_NULL */
19: PetscErrorCode PCBFBtSetOperators_BFBt(PC pc, Mat K, Mat G, Mat M)
20: {
21: PC_BFBt *ctx = (PC_BFBt *) pc->data;
24: ctx->K = K;
25: ctx->G = G;
26: ctx->M = M;
27: return(0);
28: }
33: /*@C
34: PCBFBtSetOperators - Set the matrix blocks
36: Input Parameters:
37: + pc - The PC
38: . K - The (0,0) block
39: . G - The (0,1) block
40: - M - The mass matrix associated to K, may be PETSC_NULL
42: Level: Intermediate
44: Note: K and G must not be PETSC_NULL, M can be PETSC_NULL
45: @*/
46: PetscErrorCode PCBFBtSetOperators(PC pc, Mat K, Mat G, Mat M)
47: {
48: PetscErrorCode ierr, (*f)(PC, Mat, Mat, Mat);
52: PetscObjectQueryFunction((PetscObject) pc, "PCBFBtSetOperators_C", (void (**)(void)) &f);
53: if (f) {
54: (*f)(pc, K, G, M);
55: }
56: return(0);
57: }
61: /*@C
62: PCBFBtGetKSP - Get the solver for the G^T G system
64: Input Parameter:
65: . pc - The PC
67: Output Parameter:
68: . ksp - The solver
70: Level: Intermediate
71: @*/
72: PetscErrorCode PCBFBtGetKSP(PC pc, KSP *ksp)
73: {
74: PC_BFBt *ctx = (PC_BFBt *) pc->data;
78: if (ksp) {*ksp = ctx->ksp;}
79: return(0);
80: }
84: PetscErrorCode PCBFBtCreateGtG(PC pc, Mat G, Vec inv_diag_M, Mat *GtG)
85: {
86: MPI_Comm comm;
87: Mat Ident;
88: const MatType mtype;
89: MatInfo info;
90: PetscReal fill;
91: PetscInt nnz_I, nnz_G;
92: PetscInt M, N, m, n;
96: PetscObjectGetComm((PetscObject) pc, &comm);
97: MatGetSize(G, &M, &N);
98: MatGetLocalSize(G, &m, &n);
99: MatGetType(G, &mtype);
100: MatCreate(comm, &Ident);
101: MatSetSizes(Ident, m, m, M, M);
102: MatSetType(Ident, mtype);
103: if (inv_diag_M == PETSC_NULL) {
104: Vec diag;
106: MatGetVecs(G, PETSC_NULL, &diag);
107: VecSet(diag, 1.0);
108: MatDiagonalSet(Ident, diag, INSERT_VALUES);
109: VecDestroy(diag);
110: } else {
111: MatDiagonalSet(Ident, inv_diag_M, INSERT_VALUES);
112: }
113: MatGetInfo(Ident, MAT_GLOBAL_SUM, &info);
114: nnz_I = (PetscInt) info.nz_used;
115: MatGetInfo(G, MAT_GLOBAL_SUM, &info);
116: nnz_G = (PetscInt) info.nz_used;
117: /* Not sure the best way to estimate the fill factor.
118: GtG is a laplacian on the pressure space.
119: This might tell us something useful... */
120: fill = ((PetscReal) nnz_G)/nnz_I;
121: MatPtAP(Ident, G, MAT_INITIAL_MATRIX, fill, GtG);
122: MatDestroy(Ident);
123: return(0);
124: }
126: /*
127: Performs y <- S^{-1} x
128: S^{-1} = ( G^T G )^{-1} G^T K G ( G^T G )^{-1}
129: */
132: PetscErrorCode PCApply_BFBt(PC pc, Vec x, Vec y)
133: {
134: PC_BFBt *ctx = (PC_BFBt *) pc->data;
135: KSP ksp = ctx->ksp;
136: Mat K = ctx->K;
137: Mat G = ctx->G;
138: Vec s = ctx->s;
139: Vec t = ctx->t;
140: Vec X = ctx->X;
144: /* t <- GtG_inv x */
145: KSPSolve(ksp, x, t);
146: /* s <- G t */
147: MatMult(G, t, s);
148: /* X <- K s */
149: MatMult(K, s, X);
150: /* t <- Gt X */
151: MatMultTranspose(G, X, t);
152: /* y <- GtG_inv t */
153: KSPSolve(ksp, t, y);
154: return(0);
155: }
157: /*
158: S^{-1} = ( G^T G )^{-1} G^T K G ( G^T G )^{-1}
159: = A C A
160: S^{-T} = A^T (A C)^T
161: = A^T C^T A^T, but A = G^T G which is symmetric
162: = A C^T A
163: = A G^T ( G^T K )^T A
164: = A G^T K^T G A
165: */
168: PetscErrorCode PCApplyTranspose_BFBt(PC pc, Vec x, Vec y)
169: {
170: PC_BFBt *ctx = (PC_BFBt *) pc->data;
171: KSP ksp = ctx->ksp;
172: Mat K = ctx->K;
173: Mat G = ctx->G;
174: Vec s = ctx->s;
175: Vec t = ctx->t;
176: Vec X = ctx->X;
178:
180: /* t <- GtG_inv x */
181: KSPSolve(ksp, x, t);
182: /* s <- G t */
183: MatMult(G, t, s);
184: /* X <- K^T s */
185: MatMultTranspose(K, s, X);
186: /* t <- Gt X */
187: MatMultTranspose(G, X, t);
188: /* y <- GtG_inv t */
189: KSPSolve(ksp, t, y);
190: return(0);
191: }
194: /*
195: Performs y <- S^{-1} x
196: S^{-1} = ( G^T Di G )^{-1} G^T Di K Di G ( G^T Di G )^{-1}
197: where Di = diag(M)^{-1}
198: */
201: PetscErrorCode PCApply_BFBt_diagonal_scaling(PC pc, Vec x, Vec y)
202: {
203: PC_BFBt *ctx = (PC_BFBt *) pc->data;
204: KSP ksp = ctx->ksp;
205: Mat K = ctx->K;
206: Mat G = ctx->G;
207: Vec s = ctx->s;
208: Vec t = ctx->t;
209: Vec X = ctx->X;
210: Vec di = ctx->inv_diag_M;
214: /* t <- GtG_inv x */
215: KSPSolve(ksp, x, t);
216: /* s <- G t */
217: MatMult(G, t, s);
218: /* s <- s * di */
219: VecPointwiseMult(s, s, di);
220: /* X <- K s */
221: MatMult(K, s, X);
222: /* X <- X * di */
223: VecPointwiseMult(X, X, di);
224: /* t <- Gt X */
225: MatMultTranspose(G, X, t);
226: /* y <- GtG_inv t */
227: KSPSolve(ksp, t, y);
228: return(0);
229: }
233: PetscErrorCode PCApplyTranspose_BFBt_diagonal_scaling(PC pc, Vec x, Vec y)
234: {
235: PC_BFBt *ctx = (PC_BFBt *) pc->data;
236: KSP ksp = ctx->ksp;
237: Mat K = ctx->K;
238: Mat G = ctx->G;
239: Vec s = ctx->s;
240: Vec t = ctx->t;
241: Vec X = ctx->X;
242: Vec di = ctx->inv_diag_M;
244:
246: /* t <- GtG_inv x */
247: KSPSolve(ksp, x, t);
248: /* s <- G t */
249: MatMult(G, t, s);
250: /* s <- s * di */
251: VecPointwiseMult(s, s, di);
252: /* X <- K^T s */
253: MatMultTranspose(K, s, X);
254: /* X <- X * di */
255: VecPointwiseMult(X, X, di);
256: /* t <- Gt X */
257: MatMultTranspose(G, X, t);
258: /* y <- GtG_inv t */
259: KSPSolve(ksp, t, y);
260: return(0);
261: }
265: PetscErrorCode PCSetUp_BFBt(PC pc)
266: {
267: PC_BFBt *ctx = (PC_BFBt*) pc->data;
268: MPI_Comm comm;
269: PetscTruth hasPmat;
272: PetscObjectGetComm((PetscObject) pc, &comm);
273: PCGetOperatorsSet(pc, PETSC_NULL, &hasPmat);
274: if (hasPmat) {
275: Mat pmat;
276: PetscTruth isSchur;
278: PCGetOperators(pc, PETSC_NULL, &pmat, PETSC_NULL);
279: PetscTypeCompare((PetscObject) pmat, MATSCHURCOMPLEMENT, &isSchur);
280: if (isSchur) {
281: MatSchurComplementGetSubmatrices(pmat, &ctx->K, &ctx->G, PETSC_NULL, PETSC_NULL);
282: }
283: }
284: if (ctx->K == PETSC_NULL) {SETERRQ(PETSC_ERR_SUP, "bfbt: K matrix not set");}
285: if (ctx->G == PETSC_NULL) {SETERRQ(PETSC_ERR_SUP, "bfbt: G matrix not set");}
287: /* Check for existence of objects and trash any which exist */
288: if (ctx->form_GtG && ctx->GtG != PETSC_NULL) {
289: MatDestroy(ctx->GtG);
290: ctx->GtG = PETSC_NULL;
291: }
292: if (ctx->s != PETSC_NULL) {
293: VecDestroy(ctx->s);
294: ctx->s = PETSC_NULL;
295: }
296: if (ctx->X != PETSC_NULL) {
297: VecDestroy(ctx->X);
298: ctx->X = PETSC_NULL;
299: }
300: if (ctx->t != PETSC_NULL) {
301: VecDestroy(ctx->t);
302: ctx->t = PETSC_NULL;
303: }
304: if (ctx->inv_diag_M != PETSC_NULL) {
305: VecDestroy(ctx->inv_diag_M);
306: ctx->inv_diag_M = PETSC_NULL;
307: }
308: /* Create structures */
309: MatGetVecs(ctx->K, &ctx->s, &ctx->X);
310: MatGetVecs(ctx->G, &ctx->t, PETSC_NULL);
311: if ((ctx->M != PETSC_NULL) || (ctx->scaled)) {
312: MatGetVecs(ctx->K, &ctx->inv_diag_M, PETSC_NULL);
313: if (ctx->M != PETSC_NULL) {
314: MatGetDiagonal(ctx->M, ctx->inv_diag_M);
315: VecReciprocal(ctx->inv_diag_M);
316: } else {
317: MatGetDiagonal(ctx->K, ctx->inv_diag_M);
318: VecReciprocal(ctx->inv_diag_M);
319: }
320: /* change the pc_apply routines */
321: pc->ops->apply = PCApply_BFBt_diagonal_scaling;
322: pc->ops->applytranspose = PCApplyTranspose_BFBt_diagonal_scaling;
323: }
324: PCBFBtCreateGtG(pc, ctx->G, ctx->inv_diag_M, &ctx->GtG);
325: KSPSetOperators(ctx->ksp, ctx->GtG, ctx->GtG, SAME_NONZERO_PATTERN);
326: return(0);
327: }
331: PetscErrorCode PCDestroy_BFBt(PC pc)
332: {
333: PC_BFBt *ctx = (PC_BFBt *) pc->data;
337: if (ctx == PETSC_NULL) {return(0);}
338: if (ctx->form_GtG && ctx->GtG != PETSC_NULL) {
339: MatDestroy(ctx->GtG);
340: }
341: if (ctx->ksp != PETSC_NULL) {KSPDestroy(ctx->ksp);}
342: if (ctx->s != PETSC_NULL) {VecDestroy(ctx->s);}
343: if (ctx->X != PETSC_NULL) {VecDestroy(ctx->X);}
344: if (ctx->t != PETSC_NULL) {VecDestroy(ctx->t);}
345: if (ctx->inv_diag_M != PETSC_NULL) {VecDestroy(ctx->inv_diag_M);}
346: PetscFree(ctx);
347: return(0);
348: }
352: PetscErrorCode PCSetFromOptions_BFBt(PC pc)
353: {
354: PC_BFBt *ctx = (PC_BFBt *) pc->data;
358: PetscOptionsHead("BFBt options");
359: PetscOptionsTruth("-pc_bfbt_scaled","Scale by the diagonal of K","PCBFBtSetScaled",ctx->scaled,&ctx->scaled,PETSC_NULL);
360: PetscOptionsTail();
361: KSPSetFromOptions(ctx->ksp);
362: return(0);
363: }
367: PetscErrorCode PCView_BFBt(PC pc, PetscViewer viewer)
368: {
369: PC_BFBt *ctx = (PC_BFBt *) pc->data;
373: PetscViewerASCIIPushTab(viewer);
374: if (ctx->M != PETSC_NULL) {
375: PetscViewerASCIIPrintf(viewer, "bfbt: Scaled by diag(M)\n");
376: } else if (ctx->scaled) {
377: PetscViewerASCIIPrintf(viewer, "bfbt: Scaled by diag(K)\n");
378: } else {
379: PetscViewerASCIIPrintf(viewer, "bfbt: Standard\n");
380: }
381: PetscViewerASCIIPrintf( viewer, "bfbt-ksp\n");
382: PetscViewerASCIIPushTab(viewer);
383: KSPView(ctx->ksp, viewer);
384: PetscViewerASCIIPopTab(viewer);
385: PetscViewerASCIIPopTab(viewer);
386: return(0);
387: }
392: PetscErrorCode PCCreate_BFBt(PC pc)
393: {
394: PC_BFBt *ctx = (PC_BFBt *) pc->data;
395: MPI_Comm comm;
396: const char *prefix;
400: /* create memory for ctx */
401: PetscNew(PC_BFBt, &ctx);
402: PetscLogObjectMemory(pc, sizeof(PC_BFBt));
403: /* init ctx */
404: ctx->K = PETSC_NULL;
405: ctx->G = PETSC_NULL;
406: ctx->M = PETSC_NULL;
407: ctx->inv_diag_M = PETSC_NULL;
408: ctx->GtG = PETSC_NULL;
409: ctx->s = PETSC_NULL;
410: ctx->t = PETSC_NULL;
411: ctx->X = PETSC_NULL;
412: ctx->ksp = PETSC_NULL;
413: ctx->form_GtG = PETSC_TRUE;
414: /* create internals */
415: PetscObjectGetComm((PetscObject) pc, &comm);
416: KSPCreate(comm, &ctx->ksp);
417: PCGetOptionsPrefix(pc, &prefix);
418: KSPSetOptionsPrefix(ctx->ksp, prefix);
419: KSPAppendOptionsPrefix(ctx->ksp, "pc_bfbt_");
420: /* set ctx onto pc */
421: pc->data = (void *) ctx;
422: /* define operations */
423: pc->ops->setup = PCSetUp_BFBt;
424: pc->ops->view = PCView_BFBt;
425: pc->ops->destroy = PCDestroy_BFBt;
426: pc->ops->setfromoptions = PCSetFromOptions_BFBt;
427: pc->ops->apply = PCApply_BFBt;
428: pc->ops->applytranspose = PCApplyTranspose_BFBt;
429: PetscObjectComposeFunctionDynamic((PetscObject) pc, "PCBFBtSetOperators_C", "PCBFBtSetOperators_BFBt", PCBFBtSetOperators_BFBt);
430: return(0);
431: }