Actual source code: composite.c


  2: /*
  3:       Defines a preconditioner that can consist of a collection of PCs
  4: */
  5: #include <petsc/private/pcimpl.h>
  6: #include <petscksp.h>

  8: typedef struct _PC_CompositeLink *PC_CompositeLink;
  9: struct _PC_CompositeLink {
 10:   PC               pc;
 11:   PC_CompositeLink next;
 12:   PC_CompositeLink previous;
 13: };

 15: typedef struct {
 16:   PC_CompositeLink head;
 17:   PCCompositeType  type;
 18:   Vec              work1;
 19:   Vec              work2;
 20:   PetscScalar      alpha;
 21: } PC_Composite;

 23: static PetscErrorCode PCApply_Composite_Multiplicative(PC pc, Vec x, Vec y)
 24: {
 25:   PC_Composite    *jac  = (PC_Composite *)pc->data;
 26:   PC_CompositeLink next = jac->head;
 27:   Mat              mat  = pc->pmat;



 32:   /* Set the reuse flag on children PCs */
 33:   while (next) {
 34:     PCSetReusePreconditioner(next->pc, pc->reusepreconditioner);
 35:     next = next->next;
 36:   }
 37:   next = jac->head;

 39:   if (next->next && !jac->work2) { /* allocate second work vector */
 40:     VecDuplicate(jac->work1, &jac->work2);
 41:   }
 42:   if (pc->useAmat) mat = pc->mat;
 43:   PCApply(next->pc, x, y); /* y <- B x */
 44:   while (next->next) {
 45:     next = next->next;
 46:     MatMult(mat, y, jac->work1);               /* work1 <- A y */
 47:     VecWAXPY(jac->work2, -1.0, jac->work1, x); /* work2 <- x - work1 */
 48:     PCApply(next->pc, jac->work2, jac->work1); /* work1 <- C work2 */
 49:     VecAXPY(y, 1.0, jac->work1);               /* y <- y + work1 = B x + C (x - A B x) = (B + C (1 - A B)) x */
 50:   }
 51:   if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
 52:     while (next->previous) {
 53:       next = next->previous;
 54:       MatMult(mat, y, jac->work1);
 55:       VecWAXPY(jac->work2, -1.0, jac->work1, x);
 56:       PCApply(next->pc, jac->work2, jac->work1);
 57:       VecAXPY(y, 1.0, jac->work1);
 58:     }
 59:   }
 60:   return 0;
 61: }

 63: static PetscErrorCode PCApplyTranspose_Composite_Multiplicative(PC pc, Vec x, Vec y)
 64: {
 65:   PC_Composite    *jac  = (PC_Composite *)pc->data;
 66:   PC_CompositeLink next = jac->head;
 67:   Mat              mat  = pc->pmat;

 70:   if (next->next && !jac->work2) { /* allocate second work vector */
 71:     VecDuplicate(jac->work1, &jac->work2);
 72:   }
 73:   if (pc->useAmat) mat = pc->mat;
 74:   /* locate last PC */
 75:   while (next->next) next = next->next;
 76:   PCApplyTranspose(next->pc, x, y);
 77:   while (next->previous) {
 78:     next = next->previous;
 79:     MatMultTranspose(mat, y, jac->work1);
 80:     VecWAXPY(jac->work2, -1.0, jac->work1, x);
 81:     PCApplyTranspose(next->pc, jac->work2, jac->work1);
 82:     VecAXPY(y, 1.0, jac->work1);
 83:   }
 84:   if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
 85:     next = jac->head;
 86:     while (next->next) {
 87:       next = next->next;
 88:       MatMultTranspose(mat, y, jac->work1);
 89:       VecWAXPY(jac->work2, -1.0, jac->work1, x);
 90:       PCApplyTranspose(next->pc, jac->work2, jac->work1);
 91:       VecAXPY(y, 1.0, jac->work1);
 92:     }
 93:   }
 94:   return 0;
 95: }

 97: /*
 98:     This is very special for a matrix of the form alpha I + R + S
 99: where first preconditioner is built from alpha I + S and second from
100: alpha I + R
101: */
102: static PetscErrorCode PCApply_Composite_Special(PC pc, Vec x, Vec y)
103: {
104:   PC_Composite    *jac  = (PC_Composite *)pc->data;
105:   PC_CompositeLink next = jac->head;


110:   /* Set the reuse flag on children PCs */
111:   PCSetReusePreconditioner(next->pc, pc->reusepreconditioner);
112:   PCSetReusePreconditioner(next->next->pc, pc->reusepreconditioner);

114:   PCApply(next->pc, x, jac->work1);
115:   PCApply(next->next->pc, jac->work1, y);
116:   return 0;
117: }

119: static PetscErrorCode PCApply_Composite_Additive(PC pc, Vec x, Vec y)
120: {
121:   PC_Composite    *jac  = (PC_Composite *)pc->data;
122:   PC_CompositeLink next = jac->head;


126:   /* Set the reuse flag on children PCs */
127:   while (next) {
128:     PCSetReusePreconditioner(next->pc, pc->reusepreconditioner);
129:     next = next->next;
130:   }
131:   next = jac->head;

133:   PCApply(next->pc, x, y);
134:   while (next->next) {
135:     next = next->next;
136:     PCApply(next->pc, x, jac->work1);
137:     VecAXPY(y, 1.0, jac->work1);
138:   }
139:   return 0;
140: }

142: static PetscErrorCode PCApplyTranspose_Composite_Additive(PC pc, Vec x, Vec y)
143: {
144:   PC_Composite    *jac  = (PC_Composite *)pc->data;
145:   PC_CompositeLink next = jac->head;

148:   PCApplyTranspose(next->pc, x, y);
149:   while (next->next) {
150:     next = next->next;
151:     PCApplyTranspose(next->pc, x, jac->work1);
152:     VecAXPY(y, 1.0, jac->work1);
153:   }
154:   return 0;
155: }

157: static PetscErrorCode PCSetUp_Composite(PC pc)
158: {
159:   PC_Composite    *jac  = (PC_Composite *)pc->data;
160:   PC_CompositeLink next = jac->head;
161:   DM               dm;

163:   if (!jac->work1) MatCreateVecs(pc->pmat, &jac->work1, NULL);
164:   PCGetDM(pc, &dm);
165:   while (next) {
166:     if (!next->pc->dm) PCSetDM(next->pc, dm);
167:     if (!next->pc->mat) PCSetOperators(next->pc, pc->mat, pc->pmat);
168:     next = next->next;
169:   }
170:   return 0;
171: }

173: static PetscErrorCode PCReset_Composite(PC pc)
174: {
175:   PC_Composite    *jac  = (PC_Composite *)pc->data;
176:   PC_CompositeLink next = jac->head;

178:   while (next) {
179:     PCReset(next->pc);
180:     next = next->next;
181:   }
182:   VecDestroy(&jac->work1);
183:   VecDestroy(&jac->work2);
184:   return 0;
185: }

187: static PetscErrorCode PCDestroy_Composite(PC pc)
188: {
189:   PC_Composite    *jac  = (PC_Composite *)pc->data;
190:   PC_CompositeLink next = jac->head, next_tmp;

192:   PCReset_Composite(pc);
193:   while (next) {
194:     PCDestroy(&next->pc);
195:     next_tmp = next;
196:     next     = next->next;
197:     PetscFree(next_tmp);
198:   }
199:   PetscObjectComposeFunction((PetscObject)pc, "PCCompositeSetType_C", NULL);
200:   PetscObjectComposeFunction((PetscObject)pc, "PCCompositeGetType_C", NULL);
201:   PetscObjectComposeFunction((PetscObject)pc, "PCCompositeAddPCType_C", NULL);
202:   PetscObjectComposeFunction((PetscObject)pc, "PCCompositeAddPC_C", NULL);
203:   PetscObjectComposeFunction((PetscObject)pc, "PCCompositeGetNumberPC_C", NULL);
204:   PetscObjectComposeFunction((PetscObject)pc, "PCCompositeGetPC_C", NULL);
205:   PetscObjectComposeFunction((PetscObject)pc, "PCCompositeSpecialSetAlpha_C", NULL);
206:   PetscFree(pc->data);
207:   return 0;
208: }

210: static PetscErrorCode PCSetFromOptions_Composite(PC pc, PetscOptionItems *PetscOptionsObject)
211: {
212:   PC_Composite    *jac  = (PC_Composite *)pc->data;
213:   PetscInt         nmax = 8, i;
214:   PC_CompositeLink next;
215:   char            *pcs[8];
216:   PetscBool        flg;

218:   PetscOptionsHeadBegin(PetscOptionsObject, "Composite preconditioner options");
219:   PetscOptionsEnum("-pc_composite_type", "Type of composition", "PCCompositeSetType", PCCompositeTypes, (PetscEnum)jac->type, (PetscEnum *)&jac->type, &flg);
220:   if (flg) PCCompositeSetType(pc, jac->type);
221:   PetscOptionsStringArray("-pc_composite_pcs", "List of composite solvers", "PCCompositeAddPCType", pcs, &nmax, &flg);
222:   if (flg) {
223:     for (i = 0; i < nmax; i++) {
224:       PCCompositeAddPCType(pc, pcs[i]);
225:       PetscFree(pcs[i]); /* deallocate string pcs[i], which is allocated in PetscOptionsStringArray() */
226:     }
227:   }
228:   PetscOptionsHeadEnd();

230:   next = jac->head;
231:   while (next) {
232:     PCSetFromOptions(next->pc);
233:     next = next->next;
234:   }
235:   return 0;
236: }

238: static PetscErrorCode PCView_Composite(PC pc, PetscViewer viewer)
239: {
240:   PC_Composite    *jac  = (PC_Composite *)pc->data;
241:   PC_CompositeLink next = jac->head;
242:   PetscBool        iascii;

244:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
245:   if (iascii) {
246:     PetscViewerASCIIPrintf(viewer, "Composite PC type - %s\n", PCCompositeTypes[jac->type]);
247:     PetscViewerASCIIPrintf(viewer, "PCs on composite preconditioner follow\n");
248:     PetscViewerASCIIPrintf(viewer, "---------------------------------\n");
249:   }
250:   if (iascii) PetscViewerASCIIPushTab(viewer);
251:   while (next) {
252:     PCView(next->pc, viewer);
253:     next = next->next;
254:   }
255:   if (iascii) {
256:     PetscViewerASCIIPopTab(viewer);
257:     PetscViewerASCIIPrintf(viewer, "---------------------------------\n");
258:   }
259:   return 0;
260: }

262: static PetscErrorCode PCCompositeSpecialSetAlpha_Composite(PC pc, PetscScalar alpha)
263: {
264:   PC_Composite *jac = (PC_Composite *)pc->data;

266:   jac->alpha = alpha;
267:   return 0;
268: }

270: static PetscErrorCode PCCompositeSetType_Composite(PC pc, PCCompositeType type)
271: {
272:   PC_Composite *jac = (PC_Composite *)pc->data;

274:   if (type == PC_COMPOSITE_ADDITIVE) {
275:     pc->ops->apply          = PCApply_Composite_Additive;
276:     pc->ops->applytranspose = PCApplyTranspose_Composite_Additive;
277:   } else if (type == PC_COMPOSITE_MULTIPLICATIVE || type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
278:     pc->ops->apply          = PCApply_Composite_Multiplicative;
279:     pc->ops->applytranspose = PCApplyTranspose_Composite_Multiplicative;
280:   } else if (type == PC_COMPOSITE_SPECIAL) {
281:     pc->ops->apply          = PCApply_Composite_Special;
282:     pc->ops->applytranspose = NULL;
283:   } else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Unknown composite preconditioner type");
284:   jac->type = type;
285:   return 0;
286: }

288: static PetscErrorCode PCCompositeGetType_Composite(PC pc, PCCompositeType *type)
289: {
290:   PC_Composite *jac = (PC_Composite *)pc->data;

292:   *type = jac->type;
293:   return 0;
294: }

296: static PetscErrorCode PCCompositeAddPC_Composite(PC pc, PC subpc)
297: {
298:   PC_Composite    *jac;
299:   PC_CompositeLink next, ilink;
300:   PetscInt         cnt = 0;
301:   const char      *prefix;
302:   char             newprefix[20];

304:   PetscNew(&ilink);
305:   ilink->next = NULL;
306:   ilink->pc   = subpc;

308:   jac  = (PC_Composite *)pc->data;
309:   next = jac->head;
310:   if (!next) {
311:     jac->head       = ilink;
312:     ilink->previous = NULL;
313:   } else {
314:     cnt++;
315:     while (next->next) {
316:       next = next->next;
317:       cnt++;
318:     }
319:     next->next      = ilink;
320:     ilink->previous = next;
321:   }
322:   PCGetOptionsPrefix(pc, &prefix);
323:   PCSetOptionsPrefix(subpc, prefix);
324:   PetscSNPrintf(newprefix, 20, "sub_%d_", (int)cnt);
325:   PCAppendOptionsPrefix(subpc, newprefix);
326:   PetscObjectReference((PetscObject)subpc);
327:   return 0;
328: }

330: static PetscErrorCode PCCompositeAddPCType_Composite(PC pc, PCType type)
331: {
332:   PC subpc;

334:   PCCreate(PetscObjectComm((PetscObject)pc), &subpc);
335:   PetscObjectIncrementTabLevel((PetscObject)subpc, (PetscObject)pc, 1);
336:   PCCompositeAddPC_Composite(pc, subpc);
337:   /* type is set after prefix, because some methods may modify prefix, e.g. pcksp */
338:   PCSetType(subpc, type);
339:   PCDestroy(&subpc);
340:   return 0;
341: }

343: static PetscErrorCode PCCompositeGetNumberPC_Composite(PC pc, PetscInt *n)
344: {
345:   PC_Composite    *jac;
346:   PC_CompositeLink next;

348:   jac  = (PC_Composite *)pc->data;
349:   next = jac->head;
350:   *n   = 0;
351:   while (next) {
352:     next = next->next;
353:     (*n)++;
354:   }
355:   return 0;
356: }

358: static PetscErrorCode PCCompositeGetPC_Composite(PC pc, PetscInt n, PC *subpc)
359: {
360:   PC_Composite    *jac;
361:   PC_CompositeLink next;
362:   PetscInt         i;

364:   jac  = (PC_Composite *)pc->data;
365:   next = jac->head;
366:   for (i = 0; i < n; i++) {
368:     next = next->next;
369:   }
370:   *subpc = next->pc;
371:   return 0;
372: }

374: /*@
375:    PCCompositeSetType - Sets the type of composite preconditioner.

377:    Logically Collective on pc

379:    Input Parameters:
380: +  pc - the preconditioner context
381: -  type - `PC_COMPOSITE_ADDITIVE` (default), `PC_COMPOSITE_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`

383:    Options Database Key:
384: .  -pc_composite_type <type: one of multiplicative, additive, special> - Sets composite preconditioner type

386:    Level: advanced

388: .seealso: `PCCOMPOSITE`, `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PCCompositeType`,
389:           `PCCompositeGetType()`
390: @*/
391: PetscErrorCode PCCompositeSetType(PC pc, PCCompositeType type)
392: {
395:   PetscTryMethod(pc, "PCCompositeSetType_C", (PC, PCCompositeType), (pc, type));
396:   return 0;
397: }

399: /*@
400:    PCCompositeGetType - Gets the type of composite preconditioner.

402:    Logically Collective on pc

404:    Input Parameter:
405: .  pc - the preconditioner context

407:    Output Parameter:
408: .  type - `PC_COMPOSITE_ADDITIVE` (default), `PC_COMPOSITE_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`

410:    Level: advanced

412: .seealso: `PCCOMPOSITE`, `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PCCompositeType`,
413:           `PCCompositeSetType()`
414: @*/
415: PetscErrorCode PCCompositeGetType(PC pc, PCCompositeType *type)
416: {
418:   PetscUseMethod(pc, "PCCompositeGetType_C", (PC, PCCompositeType *), (pc, type));
419:   return 0;
420: }

422: /*@
423:    PCCompositeSpecialSetAlpha - Sets alpha for the special composite preconditioner, `PC_COMPOSITE_SPECIAL`,
424:      for alphaI + R + S

426:    Logically Collective on pc

428:    Input Parameters:
429: +  pc - the preconditioner context
430: -  alpha - scale on identity

432:    Level: Developer

434: .seealso: `PCCOMPOSITE`, `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PCCompositeType`,
435:           `PCCompositeSetType()`, `PCCompositeGetType()`
436: @*/
437: PetscErrorCode PCCompositeSpecialSetAlpha(PC pc, PetscScalar alpha)
438: {
441:   PetscTryMethod(pc, "PCCompositeSpecialSetAlpha_C", (PC, PetscScalar), (pc, alpha));
442:   return 0;
443: }

445: /*@C
446:   PCCompositeAddPCType - Adds another `PC` of the given type to the composite `PC`.

448:   Collective on pc

450:   Input Parameters:
451: + pc - the preconditioner context
452: - type - the type of the new preconditioner

454:   Level: intermediate

456: .seealso: `PCCOMPOSITE`, `PCCompositeAddPC()`, `PCCompositeGetNumberPC()`
457: @*/
458: PetscErrorCode PCCompositeAddPCType(PC pc, PCType type)
459: {
461:   PetscTryMethod(pc, "PCCompositeAddPCType_C", (PC, PCType), (pc, type));
462:   return 0;
463: }

465: /*@
466:   PCCompositeAddPC - Adds another `PC` to the composite `PC`.

468:   Collective on pc

470:   Input Parameters:
471: + pc    - the preconditioner context
472: - subpc - the new preconditioner

474:    Level: intermediate

476: .seealso: `PCCOMPOSITE`, `PCCompositeAddPCType()`, `PCCompositeGetNumberPC()`
477: @*/
478: PetscErrorCode PCCompositeAddPC(PC pc, PC subpc)
479: {
482:   PetscTryMethod(pc, "PCCompositeAddPC_C", (PC, PC), (pc, subpc));
483:   return 0;
484: }

486: /*@
487:    PCCompositeGetNumberPC - Gets the number of `PC` objects in the composite `PC`.

489:    Not Collective

491:    Input Parameter:
492: .  pc - the preconditioner context

494:    Output Parameter:
495: .  num - the number of sub pcs

497:    Level: Developer

499: .seealso: `PCCOMPOSITE`, `PCCompositeGetPC()`, `PCCompositeAddPC()`,  `PCCompositeAddPCType()`
500: @*/
501: PetscErrorCode PCCompositeGetNumberPC(PC pc, PetscInt *num)
502: {
505:   PetscUseMethod(pc, "PCCompositeGetNumberPC_C", (PC, PetscInt *), (pc, num));
506:   return 0;
507: }

509: /*@
510:    PCCompositeGetPC - Gets one of the `PC` objects in the composite `PC`.

512:    Not Collective

514:    Input Parameters:
515: +  pc - the preconditioner context
516: -  n - the number of the pc requested

518:    Output Parameter:
519: .  subpc - the PC requested

521:    Level: intermediate

523:     Note:
524:     To use a different operator to construct one of the inner preconditioners first call `PCCompositeGetPC()`, then
525:     call `PCSetOperators()` on that `PC`.

527: .seealso: `PCCOMPOSITE`, `PCCompositeAddPCType()`, `PCCompositeGetNumberPC()`, `PCSetOperators()`
528: @*/
529: PetscErrorCode PCCompositeGetPC(PC pc, PetscInt n, PC *subpc)
530: {
533:   PetscUseMethod(pc, "PCCompositeGetPC_C", (PC, PetscInt, PC *), (pc, n, subpc));
534:   return 0;
535: }

537: /*MC
538:      PCCOMPOSITE - Build a preconditioner by composing together several preconditioners

540:    Options Database Keys:
541: +  -pc_composite_type <type: one of multiplicative, additive, symmetric_multiplicative, special> - Sets composite preconditioner type
542: .  -pc_use_amat - activates `PCSetUseAmat()`
543: -  -pc_composite_pcs - <pc0,pc1,...> list of PCs to compose

545:    Level: intermediate

547:    Notes:
548:    To use a Krylov method inside the composite preconditioner, set the `PCType` of one or more
549:    inner `PC`s to be `PCKSP`. Using a Krylov method inside another Krylov method can be dangerous (you get divergence or
550:    the incorrect answer) unless you use `KSPFGMRES` as the outer Krylov method

552:    To use a different operator to construct one of the inner preconditioners first call `PCCompositeGetPC()`, then
553:    call `PCSetOperators()` on that `PC`.

555: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`,
556:           `PCSHELL`, `PCKSP`, `PCCompositeSetType()`, `PCCompositeSpecialSetAlpha()`, `PCCompositeAddPCType()`,
557:           `PCCompositeGetPC()`, `PCSetUseAmat()`, `PCCompositeAddPC()`, `PCCompositeGetNumberPC()`
558: M*/

560: PETSC_EXTERN PetscErrorCode PCCreate_Composite(PC pc)
561: {
562:   PC_Composite *jac;

564:   PetscNew(&jac);

566:   pc->ops->apply           = PCApply_Composite_Additive;
567:   pc->ops->applytranspose  = PCApplyTranspose_Composite_Additive;
568:   pc->ops->setup           = PCSetUp_Composite;
569:   pc->ops->reset           = PCReset_Composite;
570:   pc->ops->destroy         = PCDestroy_Composite;
571:   pc->ops->setfromoptions  = PCSetFromOptions_Composite;
572:   pc->ops->view            = PCView_Composite;
573:   pc->ops->applyrichardson = NULL;

575:   pc->data   = (void *)jac;
576:   jac->type  = PC_COMPOSITE_ADDITIVE;
577:   jac->work1 = NULL;
578:   jac->work2 = NULL;
579:   jac->head  = NULL;

581:   PetscObjectComposeFunction((PetscObject)pc, "PCCompositeSetType_C", PCCompositeSetType_Composite);
582:   PetscObjectComposeFunction((PetscObject)pc, "PCCompositeGetType_C", PCCompositeGetType_Composite);
583:   PetscObjectComposeFunction((PetscObject)pc, "PCCompositeAddPCType_C", PCCompositeAddPCType_Composite);
584:   PetscObjectComposeFunction((PetscObject)pc, "PCCompositeAddPC_C", PCCompositeAddPC_Composite);
585:   PetscObjectComposeFunction((PetscObject)pc, "PCCompositeGetNumberPC_C", PCCompositeGetNumberPC_Composite);
586:   PetscObjectComposeFunction((PetscObject)pc, "PCCompositeGetPC_C", PCCompositeGetPC_Composite);
587:   PetscObjectComposeFunction((PetscObject)pc, "PCCompositeSpecialSetAlpha_C", PCCompositeSpecialSetAlpha_Composite);
588:   return 0;
589: }