Actual source code: snescomposite.c


  2: /*
  3:       Defines a SNES that can consist of a collection of SNESes
  4: */
  5: #include <petsc/private/snesimpl.h>
  6: #include <petscblaslapack.h>

  8: const char *const        SNESCompositeTypes[]   = {"ADDITIVE","MULTIPLICATIVE","ADDITIVEOPTIMAL","SNESCompositeType","SNES_COMPOSITE",NULL};

 10: typedef struct _SNES_CompositeLink *SNES_CompositeLink;
 11: struct _SNES_CompositeLink {
 12:   SNES               snes;
 13:   PetscReal          dmp;
 14:   Vec                X;
 15:   SNES_CompositeLink next;
 16:   SNES_CompositeLink previous;
 17: };

 19: typedef struct {
 20:   SNES_CompositeLink head;
 21:   PetscInt           nsnes;
 22:   SNESCompositeType  type;
 23:   Vec                Xorig;
 24:   PetscInt           innerFailures; /* the number of inner failures we've seen */

 26:   /* context for ADDITIVEOPTIMAL */
 27:   Vec                *Xes,*Fes;      /* solution and residual vectors for the subsolvers */
 28:   PetscReal          *fnorms;        /* norms of the solutions */
 29:   PetscScalar        *h;             /* the matrix formed as q_ij = (rdot_i, rdot_j) */
 30:   PetscScalar        *g;             /* the dotproducts of the previous function with the candidate functions */
 31:   PetscBLASInt       n;              /* matrix dimension -- nsnes */
 32:   PetscBLASInt       nrhs;           /* the number of right hand sides */
 33:   PetscBLASInt       lda;            /* the padded matrix dimension */
 34:   PetscBLASInt       ldb;            /* the padded vector dimension */
 35:   PetscReal          *s;             /* the singular values */
 36:   PetscScalar        *beta;          /* the RHS and combination */
 37:   PetscReal          rcond;          /* the exit condition */
 38:   PetscBLASInt       rank;           /* the effective rank */
 39:   PetscScalar        *work;          /* the work vector */
 40:   PetscReal          *rwork;         /* the real work vector used for complex */
 41:   PetscBLASInt       lwork;          /* the size of the work vector */
 42:   PetscBLASInt       info;           /* the output condition */

 44:   PetscReal          rtol;           /* restart tolerance for accepting the combination */
 45:   PetscReal          stol;           /* restart tolerance for the combination */
 46: } SNES_Composite;

 48: static PetscErrorCode SNESCompositeApply_Multiplicative(SNES snes,Vec X,Vec B,Vec F,PetscReal *fnorm)
 49: {
 50:   PetscErrorCode      ierr;
 51:   SNES_Composite      *jac = (SNES_Composite*)snes->data;
 52:   SNES_CompositeLink  next = jac->head;
 53:   Vec                 FSub;
 54:   SNESConvergedReason reason;

 57:   if (!next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"No composite SNESes supplied via SNESCompositeAddSNES() or -snes_composite_sneses");
 58:   if (snes->normschedule == SNES_NORM_ALWAYS) {
 59:     SNESSetInitialFunction(next->snes,F);
 60:   }
 61:   SNESSolve(next->snes,B,X);
 62:   SNESGetConvergedReason(next->snes,&reason);
 63:   if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
 64:     jac->innerFailures++;
 65:     if (jac->innerFailures >= snes->maxFailures) {
 66:       snes->reason = SNES_DIVERGED_INNER;
 67:       return(0);
 68:     }
 69:   }

 71:   while (next->next) {
 72:     /* only copy the function over in the case where the functions correspond */
 73:     if (next->snes->npcside== PC_RIGHT && next->snes->normschedule != SNES_NORM_NONE) {
 74:       SNESGetFunction(next->snes,&FSub,NULL,NULL);
 75:       next = next->next;
 76:       SNESSetInitialFunction(next->snes,FSub);
 77:     } else {
 78:       next = next->next;
 79:     }
 80:     SNESSolve(next->snes,B,X);
 81:     SNESGetConvergedReason(next->snes,&reason);
 82:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
 83:       jac->innerFailures++;
 84:       if (jac->innerFailures >= snes->maxFailures) {
 85:         snes->reason = SNES_DIVERGED_INNER;
 86:         return(0);
 87:       }
 88:     }
 89:   }
 90:   if (next->snes->npcside== PC_RIGHT) {
 91:     SNESGetFunction(next->snes,&FSub,NULL,NULL);
 92:     VecCopy(FSub,F);
 93:     if (fnorm) {
 94:       if (snes->xl && snes->xu) {
 95:         SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm);
 96:       } else {
 97:         VecNorm(F, NORM_2, fnorm);
 98:       }
 99:       SNESCheckFunctionNorm(snes,*fnorm);
100:     }
101:   } else if (snes->normschedule == SNES_NORM_ALWAYS) {
102:     SNESComputeFunction(snes,X,F);
103:     if (fnorm) {
104:       if (snes->xl && snes->xu) {
105:         SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm);
106:       } else {
107:         VecNorm(F, NORM_2, fnorm);
108:       }
109:       SNESCheckFunctionNorm(snes,*fnorm);
110:     }
111:   }
112:   return(0);
113: }

115: static PetscErrorCode SNESCompositeApply_Additive(SNES snes,Vec X,Vec B,Vec F,PetscReal *fnorm)
116: {
117:   PetscErrorCode      ierr;
118:   SNES_Composite      *jac = (SNES_Composite*)snes->data;
119:   SNES_CompositeLink  next = jac->head;
120:   Vec                 Y,Xorig;
121:   SNESConvergedReason reason;

124:   Y = snes->vec_sol_update;
125:   if (!jac->Xorig) {VecDuplicate(X,&jac->Xorig);}
126:   Xorig = jac->Xorig;
127:   VecCopy(X,Xorig);
128:   if (!next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"No composite SNESes supplied via SNESCompositeAddSNES() or -snes_composite_sneses");
129:   if (snes->normschedule == SNES_NORM_ALWAYS) {
130:     SNESSetInitialFunction(next->snes,F);
131:     while (next->next) {
132:       next = next->next;
133:       SNESSetInitialFunction(next->snes,F);
134:     }
135:   }
136:   next = jac->head;
137:   VecCopy(Xorig,Y);
138:   SNESSolve(next->snes,B,Y);
139:   SNESGetConvergedReason(next->snes,&reason);
140:   if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
141:     jac->innerFailures++;
142:     if (jac->innerFailures >= snes->maxFailures) {
143:       snes->reason = SNES_DIVERGED_INNER;
144:       return(0);
145:     }
146:   }
147:   VecAXPY(Y,-1.0,Xorig);
148:   VecAXPY(X,next->dmp,Y);
149:   while (next->next) {
150:     next = next->next;
151:     VecCopy(Xorig,Y);
152:     SNESSolve(next->snes,B,Y);
153:     SNESGetConvergedReason(next->snes,&reason);
154:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
155:       jac->innerFailures++;
156:       if (jac->innerFailures >= snes->maxFailures) {
157:         snes->reason = SNES_DIVERGED_INNER;
158:         return(0);
159:       }
160:     }
161:     VecAXPY(Y,-1.0,Xorig);
162:     VecAXPY(X,next->dmp,Y);
163:   }
164:   if (snes->normschedule == SNES_NORM_ALWAYS) {
165:     SNESComputeFunction(snes,X,F);
166:     if (fnorm) {
167:       if (snes->xl && snes->xu) {
168:         SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm);
169:       } else {
170:         VecNorm(F, NORM_2, fnorm);
171:       }
172:       SNESCheckFunctionNorm(snes,*fnorm);
173:     }
174:   }
175:   return(0);
176: }

178: /* approximately solve the overdetermined system:

180:  2*F(x_i)\cdot F(\x_j)\alpha_i = 0
181:  \alpha_i                      = 1

183:  Which minimizes the L2 norm of the linearization of:
184:  ||F(\sum_i \alpha_i*x_i)||^2

186:  With the constraint that \sum_i\alpha_i = 1
187:  Where x_i is the solution from the ith subsolver.
188:  */
189: static PetscErrorCode SNESCompositeApply_AdditiveOptimal(SNES snes,Vec X,Vec B,Vec F,PetscReal *fnorm)
190: {
191:   PetscErrorCode      ierr;
192:   SNES_Composite      *jac = (SNES_Composite*)snes->data;
193:   SNES_CompositeLink  next = jac->head;
194:   Vec                 *Xes = jac->Xes,*Fes = jac->Fes;
195:   PetscInt            i,j;
196:   PetscScalar         tot,total,ftf;
197:   PetscReal           min_fnorm;
198:   PetscInt            min_i;
199:   SNESConvergedReason reason;

202:   if (!next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"No composite SNESes supplied via SNESCompositeAddSNES() or -snes_composite_sneses");

204:   if (snes->normschedule == SNES_NORM_ALWAYS) {
205:     next = jac->head;
206:     SNESSetInitialFunction(next->snes,F);
207:     while (next->next) {
208:       next = next->next;
209:       SNESSetInitialFunction(next->snes,F);
210:     }
211:   }

213:   next = jac->head;
214:   i = 0;
215:   VecCopy(X,Xes[i]);
216:   SNESSolve(next->snes,B,Xes[i]);
217:   SNESGetConvergedReason(next->snes,&reason);
218:   if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
219:     jac->innerFailures++;
220:     if (jac->innerFailures >= snes->maxFailures) {
221:       snes->reason = SNES_DIVERGED_INNER;
222:       return(0);
223:     }
224:   }
225:   while (next->next) {
226:     i++;
227:     next = next->next;
228:     VecCopy(X,Xes[i]);
229:     SNESSolve(next->snes,B,Xes[i]);
230:     SNESGetConvergedReason(next->snes,&reason);
231:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
232:       jac->innerFailures++;
233:       if (jac->innerFailures >= snes->maxFailures) {
234:         snes->reason = SNES_DIVERGED_INNER;
235:         return(0);
236:       }
237:     }
238:   }

240:   /* all the solutions are collected; combine optimally */
241:   for (i=0;i<jac->n;i++) {
242:     for (j=0;j<i+1;j++) {
243:       VecDotBegin(Fes[i],Fes[j],&jac->h[i + j*jac->n]);
244:     }
245:     VecDotBegin(Fes[i],F,&jac->g[i]);
246:   }

248:   for (i=0;i<jac->n;i++) {
249:     for (j=0;j<i+1;j++) {
250:       VecDotEnd(Fes[i],Fes[j],&jac->h[i + j*jac->n]);
251:       if (i == j) jac->fnorms[i] = PetscSqrtReal(PetscRealPart(jac->h[i + j*jac->n]));
252:     }
253:     VecDotEnd(Fes[i],F,&jac->g[i]);
254:   }

256:   ftf = (*fnorm)*(*fnorm);

258:   for (i=0; i<jac->n; i++) {
259:     for (j=i+1;j<jac->n;j++) {
260:       jac->h[i + j*jac->n] = jac->h[j + i*jac->n];
261:     }
262:   }

264:   for (i=0; i<jac->n; i++) {
265:     for (j=0; j<jac->n; j++) {
266:       jac->h[i + j*jac->n] = jac->h[i + j*jac->n] - jac->g[j] - jac->g[i] + ftf;
267:     }
268:     jac->beta[i] = ftf - jac->g[i];
269:   }

271:   jac->info  = 0;
272:   jac->rcond = -1.;
273:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
274: #if defined(PETSC_USE_COMPLEX)
275:   PetscStackCallBLAS("LAPACKgelss",LAPACKgelss_(&jac->n,&jac->n,&jac->nrhs,jac->h,&jac->lda,jac->beta,&jac->lda,jac->s,&jac->rcond,&jac->rank,jac->work,&jac->lwork,jac->rwork,&jac->info));
276: #else
277:   PetscStackCallBLAS("LAPACKgelss",LAPACKgelss_(&jac->n,&jac->n,&jac->nrhs,jac->h,&jac->lda,jac->beta,&jac->lda,jac->s,&jac->rcond,&jac->rank,jac->work,&jac->lwork,&jac->info));
278: #endif
279:   PetscFPTrapPop();
280:   if (jac->info < 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"Bad argument to GELSS");
281:   if (jac->info > 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"SVD failed to converge");
282:   tot = 0.;
283:   total = 0.;
284:   for (i=0; i<jac->n; i++) {
285:     if (snes->errorifnotconverged && PetscIsInfOrNanScalar(jac->beta[i])) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"SVD generated inconsistent output");
286:     PetscInfo2(snes,"%D: %g\n",i,(double)PetscRealPart(jac->beta[i]));
287:     tot += jac->beta[i];
288:     total += PetscAbsScalar(jac->beta[i]);
289:   }
290:   VecScale(X,(1. - tot));
291:   VecMAXPY(X,jac->n,jac->beta,Xes);
292:   SNESComputeFunction(snes,X,F);

294:   if (snes->xl && snes->xu) {
295:     SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm);
296:   } else {
297:     VecNorm(F, NORM_2, fnorm);
298:   }

300:   /* take the minimum-normed candidate if it beats the combination by a factor of rtol or the combination has stagnated */
301:   min_fnorm = jac->fnorms[0];
302:   min_i     = 0;
303:   for (i=0; i<jac->n; i++) {
304:     if (jac->fnorms[i] < min_fnorm) {
305:       min_fnorm = jac->fnorms[i];
306:       min_i     = i;
307:     }
308:   }

310:   /* stagnation or divergence restart to the solution of the solver that failed the least */
311:   if (PetscRealPart(total) < jac->stol || min_fnorm*jac->rtol < *fnorm) {
312:     VecCopy(jac->Xes[min_i],X);
313:     VecCopy(jac->Fes[min_i],F);
314:     *fnorm = min_fnorm;
315:   }
316:   return(0);
317: }

319: static PetscErrorCode SNESSetUp_Composite(SNES snes)
320: {
321:   PetscErrorCode     ierr;
322:   DM                 dm;
323:   SNES_Composite     *jac = (SNES_Composite*)snes->data;
324:   SNES_CompositeLink next = jac->head;
325:   PetscInt           n=0,i;
326:   Vec                F;

329:   SNESGetDM(snes,&dm);

331:   if (snes->ops->computevariablebounds) {
332:     /* SNESVI only ever calls computevariablebounds once, so calling it once here is justified */
333:     if (!snes->xl) {VecDuplicate(snes->vec_sol,&snes->xl);}
334:     if (!snes->xu) {VecDuplicate(snes->vec_sol,&snes->xu);}
335:     (*snes->ops->computevariablebounds)(snes,snes->xl,snes->xu);
336:   }

338:   while (next) {
339:     n++;
340:     SNESSetDM(next->snes,dm);
341:     SNESSetApplicationContext(next->snes, snes->user);
342:     if (snes->xl && snes->xu) {
343:       if (snes->ops->computevariablebounds) {
344:         SNESVISetComputeVariableBounds(next->snes, snes->ops->computevariablebounds);
345:       } else {
346:         SNESVISetVariableBounds(next->snes,snes->xl,snes->xu);
347:       }
348:     }

350:     next = next->next;
351:   }
352:   jac->nsnes = n;
353:   SNESGetFunction(snes,&F,NULL,NULL);
354:   if (jac->type == SNES_COMPOSITE_ADDITIVEOPTIMAL) {
355:     VecDuplicateVecs(F,jac->nsnes,&jac->Xes);
356:     PetscMalloc1(n,&jac->Fes);
357:     PetscMalloc1(n,&jac->fnorms);
358:     next = jac->head;
359:     i = 0;
360:     while (next) {
361:       SNESGetFunction(next->snes,&F,NULL,NULL);
362:       jac->Fes[i] = F;
363:       PetscObjectReference((PetscObject)F);
364:       next = next->next;
365:       i++;
366:     }
367:     /* allocate the subspace direct solve area */
368:     jac->nrhs  = 1;
369:     jac->lda   = jac->nsnes;
370:     jac->ldb   = jac->nsnes;
371:     jac->n     = jac->nsnes;

373:     PetscMalloc1(jac->n*jac->n,&jac->h);
374:     PetscMalloc1(jac->n,&jac->beta);
375:     PetscMalloc1(jac->n,&jac->s);
376:     PetscMalloc1(jac->n,&jac->g);
377:     jac->lwork = 12*jac->n;
378: #if defined(PETSC_USE_COMPLEX)
379:     PetscMalloc1(jac->lwork,&jac->rwork);
380: #endif
381:     PetscMalloc1(jac->lwork,&jac->work);
382:   }

384:   return(0);
385: }

387: static PetscErrorCode SNESReset_Composite(SNES snes)
388: {
389:   SNES_Composite     *jac = (SNES_Composite*)snes->data;
390:   PetscErrorCode   ierr;
391:   SNES_CompositeLink next = jac->head;

394:   while (next) {
395:     SNESReset(next->snes);
396:     next = next->next;
397:   }
398:   VecDestroy(&jac->Xorig);
399:   if (jac->Xes) {VecDestroyVecs(jac->nsnes,&jac->Xes);}
400:   if (jac->Fes) {VecDestroyVecs(jac->nsnes,&jac->Fes);}
401:   PetscFree(jac->fnorms);
402:   PetscFree(jac->h);
403:   PetscFree(jac->s);
404:   PetscFree(jac->g);
405:   PetscFree(jac->beta);
406:   PetscFree(jac->work);
407:   PetscFree(jac->rwork);
408:   return(0);
409: }

411: static PetscErrorCode SNESDestroy_Composite(SNES snes)
412: {
413:   SNES_Composite     *jac = (SNES_Composite*)snes->data;
414:   PetscErrorCode     ierr;
415:   SNES_CompositeLink next = jac->head,next_tmp;

418:   SNESReset_Composite(snes);
419:   while (next) {
420:     SNESDestroy(&next->snes);
421:     next_tmp = next;
422:     next     = next->next;
423:     PetscFree(next_tmp);
424:   }
425:   PetscFree(snes->data);
426:   return(0);
427: }

429: static PetscErrorCode SNESSetFromOptions_Composite(PetscOptionItems *PetscOptionsObject,SNES snes)
430: {
431:   SNES_Composite     *jac = (SNES_Composite*)snes->data;
432:   PetscErrorCode     ierr;
433:   PetscInt           nmax = 8,i;
434:   SNES_CompositeLink next;
435:   char               *sneses[8];
436:   PetscReal          dmps[8];
437:   PetscBool          flg;

440:   PetscOptionsHead(PetscOptionsObject,"Composite preconditioner options");
441:   PetscOptionsEnum("-snes_composite_type","Type of composition","SNESCompositeSetType",SNESCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&jac->type,&flg);
442:   if (flg) {
443:     SNESCompositeSetType(snes,jac->type);
444:   }
445:   PetscOptionsStringArray("-snes_composite_sneses","List of composite solvers","SNESCompositeAddSNES",sneses,&nmax,&flg);
446:   if (flg) {
447:     for (i=0; i<nmax; i++) {
448:       SNESCompositeAddSNES(snes,sneses[i]);
449:       PetscFree(sneses[i]);   /* deallocate string sneses[i], which is allocated in PetscOptionsStringArray() */
450:     }
451:   }
452:   PetscOptionsRealArray("-snes_composite_damping","Damping of the additive composite solvers","SNESCompositeSetDamping",dmps,&nmax,&flg);
453:   if (flg) {
454:     for (i=0; i<nmax; i++) {
455:       SNESCompositeSetDamping(snes,i,dmps[i]);
456:     }
457:   }
458:   PetscOptionsReal("-snes_composite_stol","Step tolerance for restart on the additive composite solvers","",jac->stol,&jac->stol,NULL);
459:   PetscOptionsReal("-snes_composite_rtol","Residual tolerance for the additive composite solvers","",jac->rtol,&jac->rtol,NULL);
460:   PetscOptionsTail();

462:   next = jac->head;
463:   while (next) {
464:     SNESSetFromOptions(next->snes);
465:     next = next->next;
466:   }
467:   return(0);
468: }

470: static PetscErrorCode SNESView_Composite(SNES snes,PetscViewer viewer)
471: {
472:   SNES_Composite     *jac = (SNES_Composite*)snes->data;
473:   PetscErrorCode     ierr;
474:   SNES_CompositeLink next = jac->head;
475:   PetscBool          iascii;

478:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
479:   if (iascii) {
480:     PetscViewerASCIIPrintf(viewer,"  type - %s\n",SNESCompositeTypes[jac->type]);
481:     PetscViewerASCIIPrintf(viewer,"  SNESes on composite preconditioner follow\n");
482:     PetscViewerASCIIPrintf(viewer,"  ---------------------------------\n");
483:   }
484:   if (iascii) {
485:     PetscViewerASCIIPushTab(viewer);
486:   }
487:   while (next) {
488:     SNESView(next->snes,viewer);
489:     next = next->next;
490:   }
491:   if (iascii) {
492:     PetscViewerASCIIPopTab(viewer);
493:     PetscViewerASCIIPrintf(viewer,"  ---------------------------------\n");
494:   }
495:   return(0);
496: }

498: /* ------------------------------------------------------------------------------*/

500: static PetscErrorCode  SNESCompositeSetType_Composite(SNES snes,SNESCompositeType type)
501: {
502:   SNES_Composite *jac = (SNES_Composite*)snes->data;

505:   jac->type = type;
506:   return(0);
507: }

509: static PetscErrorCode  SNESCompositeAddSNES_Composite(SNES snes,SNESType type)
510: {
511:   SNES_Composite     *jac;
512:   SNES_CompositeLink next,ilink;
513:   PetscErrorCode     ierr;
514:   PetscInt           cnt = 0;
515:   const char         *prefix;
516:   char               newprefix[20];
517:   DM                 dm;

520:   PetscNewLog(snes,&ilink);
521:   ilink->next = NULL;
522:   SNESCreate(PetscObjectComm((PetscObject)snes),&ilink->snes);
523:   PetscLogObjectParent((PetscObject)snes,(PetscObject)ilink->snes);
524:   SNESGetDM(snes,&dm);
525:   SNESSetDM(ilink->snes,dm);
526:   SNESSetTolerances(ilink->snes,snes->abstol,snes->rtol,snes->stol,1,snes->max_funcs);
527:   PetscObjectCopyFortranFunctionPointers((PetscObject)snes,(PetscObject)ilink->snes);
528:   jac  = (SNES_Composite*)snes->data;
529:   next = jac->head;
530:   if (!next) {
531:     jac->head       = ilink;
532:     ilink->previous = NULL;
533:   } else {
534:     cnt++;
535:     while (next->next) {
536:       next = next->next;
537:       cnt++;
538:     }
539:     next->next      = ilink;
540:     ilink->previous = next;
541:   }
542:   SNESGetOptionsPrefix(snes,&prefix);
543:   SNESSetOptionsPrefix(ilink->snes,prefix);
544:   PetscSNPrintf(newprefix,sizeof(newprefix),"sub_%d_",(int)cnt);
545:   SNESAppendOptionsPrefix(ilink->snes,newprefix);
546:   PetscObjectIncrementTabLevel((PetscObject)ilink->snes,(PetscObject)snes,1);
547:   SNESSetType(ilink->snes,type);
548:   SNESSetNormSchedule(ilink->snes, SNES_NORM_FINAL_ONLY);

550:   ilink->dmp = 1.0;
551:   jac->nsnes++;
552:   return(0);
553: }

555: static PetscErrorCode  SNESCompositeGetSNES_Composite(SNES snes,PetscInt n,SNES *subsnes)
556: {
557:   SNES_Composite     *jac;
558:   SNES_CompositeLink next;
559:   PetscInt           i;

562:   jac  = (SNES_Composite*)snes->data;
563:   next = jac->head;
564:   for (i=0; i<n; i++) {
565:     if (!next->next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Not enough SNESes in composite preconditioner");
566:     next = next->next;
567:   }
568:   *subsnes = next->snes;
569:   return(0);
570: }

572: /* -------------------------------------------------------------------------------- */
573: /*@C
574:    SNESCompositeSetType - Sets the type of composite preconditioner.

576:    Logically Collective on SNES

578:    Input Parameters:
579: +  snes - the preconditioner context
580: -  type - SNES_COMPOSITE_ADDITIVE (default), SNES_COMPOSITE_MULTIPLICATIVE

582:    Options Database Key:
583: .  -snes_composite_type <type: one of multiplicative, additive, special> - Sets composite preconditioner type

585:    Level: Developer

587: @*/
588: PetscErrorCode  SNESCompositeSetType(SNES snes,SNESCompositeType type)
589: {

595:   PetscTryMethod(snes,"SNESCompositeSetType_C",(SNES,SNESCompositeType),(snes,type));
596:   return(0);
597: }

599: /*@C
600:    SNESCompositeAddSNES - Adds another SNES to the composite SNES.

602:    Collective on SNES

604:    Input Parameters:
605: +  snes - the preconditioner context
606: -  type - the type of the new preconditioner

608:    Level: Developer

610: @*/
611: PetscErrorCode  SNESCompositeAddSNES(SNES snes,SNESType type)
612: {

617:   PetscTryMethod(snes,"SNESCompositeAddSNES_C",(SNES,SNESType),(snes,type));
618:   return(0);
619: }

621: /*@
622:    SNESCompositeGetSNES - Gets one of the SNES objects in the composite SNES.

624:    Not Collective

626:    Input Parameters:
627: +  snes - the preconditioner context
628: -  n - the number of the snes requested

630:    Output Parameters:
631: .  subsnes - the SNES requested

633:    Level: Developer

635: .seealso: SNESCompositeAddSNES()
636: @*/
637: PetscErrorCode  SNESCompositeGetSNES(SNES snes,PetscInt n,SNES *subsnes)
638: {

644:   PetscUseMethod(snes,"SNESCompositeGetSNES_C",(SNES,PetscInt,SNES*),(snes,n,subsnes));
645:   return(0);
646: }

648: /*@
649:    SNESCompositeGetNumber - Get the number of subsolvers in the composite SNES.

651:    Logically Collective on SNES

653:    Input Parameter:
654:    snes - the preconditioner context

656:    Output Parameter:
657:    n - the number of subsolvers

659:    Level: Developer

661: @*/
662: PetscErrorCode  SNESCompositeGetNumber(SNES snes,PetscInt *n)
663: {
664:   SNES_Composite     *jac;
665:   SNES_CompositeLink next;

668:   jac  = (SNES_Composite*)snes->data;
669:   next = jac->head;

671:   *n = 0;
672:   while (next) {
673:     *n = *n + 1;
674:     next = next->next;
675:   }
676:   return(0);
677: }

679: static PetscErrorCode  SNESCompositeSetDamping_Composite(SNES snes,PetscInt n,PetscReal dmp)
680: {
681:   SNES_Composite     *jac;
682:   SNES_CompositeLink next;
683:   PetscInt           i;

686:   jac  = (SNES_Composite*)snes->data;
687:   next = jac->head;
688:   for (i=0; i<n; i++) {
689:     if (!next->next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Not enough SNESes in composite preconditioner");
690:     next = next->next;
691:   }
692:   next->dmp = dmp;
693:   return(0);
694: }

696: /*@
697:    SNESCompositeSetDamping - Sets the damping of a subsolver when using additive composite SNES.

699:    Not Collective

701:    Input Parameters:
702: +  snes - the preconditioner context
703: .  n - the number of the snes requested
704: -  dmp - the damping

706:    Level: Developer

708: .seealso: SNESCompositeAddSNES()
709: @*/
710: PetscErrorCode  SNESCompositeSetDamping(SNES snes,PetscInt n,PetscReal dmp)
711: {

716:   PetscUseMethod(snes,"SNESCompositeSetDamping_C",(SNES,PetscInt,PetscReal),(snes,n,dmp));
717:   return(0);
718: }

720: static PetscErrorCode SNESSolve_Composite(SNES snes)
721: {
722:   Vec              F,X,B,Y;
723:   PetscInt         i;
724:   PetscReal        fnorm = 0.0, xnorm = 0.0, snorm = 0.0;
725:   PetscErrorCode   ierr;
726:   SNESNormSchedule normtype;
727:   SNES_Composite   *comp = (SNES_Composite*)snes->data;

730:   X = snes->vec_sol;
731:   F = snes->vec_func;
732:   B = snes->vec_rhs;
733:   Y = snes->vec_sol_update;

735:   PetscObjectSAWsTakeAccess((PetscObject)snes);
736:   snes->iter   = 0;
737:   snes->norm   = 0.;
738:   comp->innerFailures = 0;
739:   PetscObjectSAWsGrantAccess((PetscObject)snes);
740:   snes->reason = SNES_CONVERGED_ITERATING;
741:   SNESGetNormSchedule(snes, &normtype);
742:   if (normtype == SNES_NORM_ALWAYS || normtype == SNES_NORM_INITIAL_ONLY || normtype == SNES_NORM_INITIAL_FINAL_ONLY) {
743:     if (!snes->vec_func_init_set) {
744:       SNESComputeFunction(snes,X,F);
745:     } else snes->vec_func_init_set = PETSC_FALSE;

747:     if (snes->xl && snes->xu) {
748:       SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm);
749:     } else {
750:       VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
751:     }
752:     SNESCheckFunctionNorm(snes,fnorm);
753:     PetscObjectSAWsTakeAccess((PetscObject)snes);
754:     snes->iter = 0;
755:     snes->norm = fnorm;
756:     PetscObjectSAWsGrantAccess((PetscObject)snes);
757:     SNESLogConvergenceHistory(snes,snes->norm,0);
758:     SNESMonitor(snes,0,snes->norm);

760:     /* test convergence */
761:     (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
762:     if (snes->reason) return(0);
763:   } else {
764:     PetscObjectSAWsGrantAccess((PetscObject)snes);
765:     SNESLogConvergenceHistory(snes,snes->norm,0);
766:     SNESMonitor(snes,0,snes->norm);
767:   }

769:   for (i = 0; i < snes->max_its; i++) {
770:     /* Call general purpose update function */
771:     if (snes->ops->update) {
772:       (*snes->ops->update)(snes, snes->iter);
773:     }

775:     /* Copy the state before modification by application of the composite solver;
776:        we will subtract the new state after application */
777:     VecCopy(X, Y);

779:     if (comp->type == SNES_COMPOSITE_ADDITIVE) {
780:       SNESCompositeApply_Additive(snes,X,B,F,&fnorm);
781:     } else if (comp->type == SNES_COMPOSITE_MULTIPLICATIVE) {
782:       SNESCompositeApply_Multiplicative(snes,X,B,F,&fnorm);
783:     } else if (comp->type == SNES_COMPOSITE_ADDITIVEOPTIMAL) {
784:       SNESCompositeApply_AdditiveOptimal(snes,X,B,F,&fnorm);
785:     } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Unsupported SNESComposite type");
786:     if (snes->reason < 0) break;

788:     /* Compute the solution update for convergence testing */
789:     VecAYPX(Y, -1.0, X);

791:     if ((i == snes->max_its - 1) && (normtype == SNES_NORM_INITIAL_FINAL_ONLY || normtype == SNES_NORM_FINAL_ONLY)) {
792:       SNESComputeFunction(snes,X,F);

794:       if (snes->xl && snes->xu) {
795:         VecNormBegin(X, NORM_2, &xnorm);
796:         VecNormBegin(Y, NORM_2, &snorm);
797:         SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm);
798:         VecNormEnd(X, NORM_2, &xnorm);
799:         VecNormEnd(Y, NORM_2, &snorm);
800:       } else {
801:         VecNormBegin(F, NORM_2, &fnorm);
802:         VecNormBegin(X, NORM_2, &xnorm);
803:         VecNormBegin(Y, NORM_2, &snorm);

805:         VecNormEnd(F, NORM_2, &fnorm);
806:         VecNormEnd(X, NORM_2, &xnorm);
807:         VecNormEnd(Y, NORM_2, &snorm);
808:       }
809:       SNESCheckFunctionNorm(snes,fnorm);
810:     } else if (normtype == SNES_NORM_ALWAYS) {
811:       VecNormBegin(X, NORM_2, &xnorm);
812:       VecNormBegin(Y, NORM_2, &snorm);
813:       VecNormEnd(X, NORM_2, &xnorm);
814:       VecNormEnd(Y, NORM_2, &snorm);
815:     }
816:     /* Monitor convergence */
817:     PetscObjectSAWsTakeAccess((PetscObject)snes);
818:     snes->iter = i+1;
819:     snes->norm = fnorm;
820:     snes->xnorm = xnorm;
821:     snes->ynorm = snorm;
822:     PetscObjectSAWsGrantAccess((PetscObject)snes);
823:     SNESLogConvergenceHistory(snes,snes->norm,0);
824:     SNESMonitor(snes,snes->iter,snes->norm);
825:     /* Test for convergence */
826:     if (normtype == SNES_NORM_ALWAYS) {(*snes->ops->converged)(snes,snes->iter,xnorm,snorm,fnorm,&snes->reason,snes->cnvP);}
827:     if (snes->reason) break;
828:   }
829:   if (normtype == SNES_NORM_ALWAYS) {
830:     if (i == snes->max_its) {
831:       PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);
832:       if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
833:     }
834:   } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS;
835:   return(0);
836: }

838: /* -------------------------------------------------------------------------------------------*/

840: /*MC
841:      SNESCOMPOSITE - Build a preconditioner by composing together several nonlinear solvers

843:    Options Database Keys:
844: +  -snes_composite_type <type: one of multiplicative, additive, symmetric_multiplicative, special> - Sets composite preconditioner type
845: -  -snes_composite_sneses - <snes0,snes1,...> list of SNESes to compose

847:    Level: intermediate

849: .seealso:  SNESCreate(), SNESSetType(), SNESType (for list of available types), SNES,
850:            SNESSHELL, SNESCompositeSetType(), SNESCompositeSpecialSetAlpha(), SNESCompositeAddSNES(),
851:            SNESCompositeGetSNES()

853:    References:
854: .  1. - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
855:    SIAM Review, 57(4), 2015

857: M*/

859: PETSC_EXTERN PetscErrorCode SNESCreate_Composite(SNES snes)
860: {
862:   SNES_Composite *jac;

865:   PetscNewLog(snes,&jac);

867:   snes->ops->solve           = SNESSolve_Composite;
868:   snes->ops->setup           = SNESSetUp_Composite;
869:   snes->ops->reset           = SNESReset_Composite;
870:   snes->ops->destroy         = SNESDestroy_Composite;
871:   snes->ops->setfromoptions  = SNESSetFromOptions_Composite;
872:   snes->ops->view            = SNESView_Composite;

874:   snes->usesksp        = PETSC_FALSE;

876:   snes->alwayscomputesfinalresidual = PETSC_FALSE;

878:   snes->data = (void*)jac;
879:   jac->type  = SNES_COMPOSITE_ADDITIVEOPTIMAL;
880:   jac->Fes   = NULL;
881:   jac->Xes   = NULL;
882:   jac->fnorms = NULL;
883:   jac->nsnes = 0;
884:   jac->head  = NULL;
885:   jac->stol  = 0.1;
886:   jac->rtol  = 1.1;

888:   jac->h     = NULL;
889:   jac->s     = NULL;
890:   jac->beta  = NULL;
891:   jac->work  = NULL;
892:   jac->rwork = NULL;

894:   PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeSetType_C",SNESCompositeSetType_Composite);
895:   PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeAddSNES_C",SNESCompositeAddSNES_Composite);
896:   PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeGetSNES_C",SNESCompositeGetSNES_Composite);
897:   PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeSetDamping_C",SNESCompositeSetDamping_Composite);
898:   return(0);
899: }