Actual source code: snescomposite.c

petsc-3.9.4 2018-09-11
Report Typos and Errors

  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",0};

 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: #if defined(PETSC_MISSING_LAPACK_GELSS)
272:   SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"SNESCOMPOSITE with ADDITIVEOPTIMAL requires the LAPACK GELSS routine.");
273: #else
274:   jac->info  = 0;
275:   jac->rcond = -1.;
276:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
277: #if defined(PETSC_USE_COMPLEX)
278:   PetscStackCall("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));
279: #else
280:   PetscStackCall("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));
281: #endif
282:   PetscFPTrapPop();
283:   if (jac->info < 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"Bad argument to GELSS");
284:   if (jac->info > 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"SVD failed to converge");
285: #endif
286:   tot = 0.;
287:   total = 0.;
288:   for (i=0; i<jac->n; i++) {
289:     if (snes->errorifnotconverged && PetscIsInfOrNanScalar(jac->beta[i])) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"SVD generated inconsistent output");
290:     PetscInfo2(snes,"%D: %g\n",i,(double)PetscRealPart(jac->beta[i]));
291:     tot += jac->beta[i];
292:     total += PetscAbsScalar(jac->beta[i]);
293:   }
294:   VecScale(X,(1. - tot));
295:   VecMAXPY(X,jac->n,jac->beta,Xes);
296:   SNESComputeFunction(snes,X,F);

298:   if (snes->xl && snes->xu) {
299:     SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm);
300:   } else {
301:     VecNorm(F, NORM_2, fnorm);
302:   }

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

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

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

333:   SNESGetDM(snes,&dm);

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

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

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

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

388:   return(0);
389: }

391: static PetscErrorCode SNESReset_Composite(SNES snes)
392: {
393:   SNES_Composite     *jac = (SNES_Composite*)snes->data;
394:   PetscErrorCode   ierr;
395:   SNES_CompositeLink next = jac->head;

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

415: static PetscErrorCode SNESDestroy_Composite(SNES snes)
416: {
417:   SNES_Composite     *jac = (SNES_Composite*)snes->data;
418:   PetscErrorCode     ierr;
419:   SNES_CompositeLink next = jac->head,next_tmp;

422:   SNESReset_Composite(snes);
423:   while (next) {
424:     SNESDestroy(&next->snes);
425:     next_tmp = next;
426:     next     = next->next;
427:     PetscFree(next_tmp);
428:   }
429:   PetscFree(snes->data);
430:   return(0);
431: }

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

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

466:   next = jac->head;
467:   while (next) {
468:     SNESSetFromOptions(next->snes);
469:     next = next->next;
470:   }
471:   return(0);
472: }

474: static PetscErrorCode SNESView_Composite(SNES snes,PetscViewer viewer)
475: {
476:   SNES_Composite     *jac = (SNES_Composite*)snes->data;
477:   PetscErrorCode     ierr;
478:   SNES_CompositeLink next = jac->head;
479:   PetscBool          iascii;

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

502: /* ------------------------------------------------------------------------------*/

504: static PetscErrorCode  SNESCompositeSetType_Composite(SNES snes,SNESCompositeType type)
505: {
506:   SNES_Composite *jac = (SNES_Composite*)snes->data;

509:   jac->type = type;
510:   return(0);
511: }

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

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

554:   ilink->dmp = 1.0;
555:   jac->nsnes++;
556:   return(0);
557: }

559: static PetscErrorCode  SNESCompositeGetSNES_Composite(SNES snes,PetscInt n,SNES *subsnes)
560: {
561:   SNES_Composite     *jac;
562:   SNES_CompositeLink next;
563:   PetscInt           i;

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

576: /* -------------------------------------------------------------------------------- */
577: /*@C
578:    SNESCompositeSetType - Sets the type of composite preconditioner.

580:    Logically Collective on SNES

582:    Input Parameter:
583: +  snes - the preconditioner context
584: -  type - SNES_COMPOSITE_ADDITIVE (default), SNES_COMPOSITE_MULTIPLICATIVE

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

589:    Level: Developer

591: .keywords: SNES, set, type, composite preconditioner, additive, multiplicative
592: @*/
593: PetscErrorCode  SNESCompositeSetType(SNES snes,SNESCompositeType type)
594: {

600:   PetscTryMethod(snes,"SNESCompositeSetType_C",(SNES,SNESCompositeType),(snes,type));
601:   return(0);
602: }

604: /*@C
605:    SNESCompositeAddSNES - Adds another SNES to the composite SNES.

607:    Collective on SNES

609:    Input Parameters:
610: +  snes - the preconditioner context
611: -  type - the type of the new preconditioner

613:    Level: Developer

615: .keywords: SNES, composite preconditioner, add
616: @*/
617: PetscErrorCode  SNESCompositeAddSNES(SNES snes,SNESType type)
618: {

623:   PetscTryMethod(snes,"SNESCompositeAddSNES_C",(SNES,SNESType),(snes,type));
624:   return(0);
625: }
626: /*@
627:    SNESCompositeGetSNES - Gets one of the SNES objects in the composite SNES.

629:    Not Collective

631:    Input Parameter:
632: +  snes - the preconditioner context
633: -  n - the number of the snes requested

635:    Output Parameters:
636: .  subsnes - the SNES requested

638:    Level: Developer

640: .keywords: SNES, get, composite preconditioner, sub preconditioner

642: .seealso: SNESCompositeAddSNES()
643: @*/
644: PetscErrorCode  SNESCompositeGetSNES(SNES snes,PetscInt n,SNES *subsnes)
645: {

651:   PetscUseMethod(snes,"SNESCompositeGetSNES_C",(SNES,PetscInt,SNES*),(snes,n,subsnes));
652:   return(0);
653: }

655: /*@
656:    SNESCompositeGetNumber - Get the number of subsolvers in the composite SNES.

658:    Logically Collective on SNES

660:    Input Parameter:
661:    snes - the preconditioner context

663:    Output Parameter:
664:    n - the number of subsolvers

666:    Level: Developer

668: .keywords: SNES, composite preconditioner
669: @*/
670: PetscErrorCode  SNESCompositeGetNumber(SNES snes,PetscInt *n)
671: {
672:   SNES_Composite     *jac;
673:   SNES_CompositeLink next;

676:   jac  = (SNES_Composite*)snes->data;
677:   next = jac->head;

679:   *n = 0;
680:   while (next) {
681:     *n = *n + 1;
682:     next = next->next;
683:   }
684:   return(0);
685: }

687: static PetscErrorCode  SNESCompositeSetDamping_Composite(SNES snes,PetscInt n,PetscReal dmp)
688: {
689:   SNES_Composite     *jac;
690:   SNES_CompositeLink next;
691:   PetscInt           i;

694:   jac  = (SNES_Composite*)snes->data;
695:   next = jac->head;
696:   for (i=0; i<n; i++) {
697:     if (!next->next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Not enough SNESes in composite preconditioner");
698:     next = next->next;
699:   }
700:   next->dmp = dmp;
701:   return(0);
702: }

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

707:    Not Collective

709:    Input Parameter:
710: +  snes - the preconditioner context
711: .  n - the number of the snes requested
712: -  dmp - the damping

714:    Level: Developer

716: .keywords: SNES, get, composite preconditioner, sub preconditioner

718: .seealso: SNESCompositeAddSNES()
719: @*/
720: PetscErrorCode  SNESCompositeSetDamping(SNES snes,PetscInt n,PetscReal dmp)
721: {

726:   PetscUseMethod(snes,"SNESCompositeSetDamping_C",(SNES,PetscInt,PetscReal),(snes,n,dmp));
727:   return(0);
728: }

730: static PetscErrorCode SNESSolve_Composite(SNES snes)
731: {
732:   Vec              F;
733:   Vec              X;
734:   Vec              B;
735:   PetscInt         i;
736:   PetscReal        fnorm = 0.0, xnorm = 0.0, snorm = 0.0;
737:   PetscErrorCode   ierr;
738:   SNESNormSchedule normtype;
739:   SNES_Composite   *comp = (SNES_Composite*)snes->data;

742:   X = snes->vec_sol;
743:   F = snes->vec_func;
744:   B = snes->vec_rhs;

746:   PetscObjectSAWsTakeAccess((PetscObject)snes);
747:   snes->iter   = 0;
748:   snes->norm   = 0.;
749:   comp->innerFailures = 0;
750:   PetscObjectSAWsGrantAccess((PetscObject)snes);
751:   SNESSetWorkVecs(snes, 1);
752:   snes->reason = SNES_CONVERGED_ITERATING;
753:   SNESGetNormSchedule(snes, &normtype);
754:   if (normtype == SNES_NORM_ALWAYS || normtype == SNES_NORM_INITIAL_ONLY || normtype == SNES_NORM_INITIAL_FINAL_ONLY) {
755:     if (!snes->vec_func_init_set) {
756:       SNESComputeFunction(snes,X,F);
757:     } else snes->vec_func_init_set = PETSC_FALSE;

759:     if (snes->xl && snes->xu) {
760:       SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm);
761:     } else {
762:       VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
763:     }
764:     SNESCheckFunctionNorm(snes,fnorm);
765:     PetscObjectSAWsTakeAccess((PetscObject)snes);
766:     snes->iter = 0;
767:     snes->norm = fnorm;
768:     PetscObjectSAWsGrantAccess((PetscObject)snes);
769:     SNESLogConvergenceHistory(snes,snes->norm,0);
770:     SNESMonitor(snes,0,snes->norm);

772:     /* test convergence */
773:     (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
774:     if (snes->reason) return(0);
775:   } else {
776:     PetscObjectSAWsGrantAccess((PetscObject)snes);
777:     SNESLogConvergenceHistory(snes,snes->norm,0);
778:     SNESMonitor(snes,0,snes->norm);
779:   }

781:   /* Call general purpose update function */
782:   if (snes->ops->update) {
783:     (*snes->ops->update)(snes, snes->iter);
784:   }

786:   for (i = 0; i < snes->max_its; i++) {
787:     /* Copy the state before modification by application of the composite solver;
788:        we will subtract the new state after application */
789:     VecCopy(X, snes->work[0]);

791:     if (comp->type == SNES_COMPOSITE_ADDITIVE) {
792:       SNESCompositeApply_Additive(snes,X,B,F,&fnorm);
793:     } else if (comp->type == SNES_COMPOSITE_MULTIPLICATIVE) {
794:       SNESCompositeApply_Multiplicative(snes,X,B,F,&fnorm);
795:     } else if (comp->type == SNES_COMPOSITE_ADDITIVEOPTIMAL) {
796:       SNESCompositeApply_AdditiveOptimal(snes,X,B,F,&fnorm);
797:     } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Unsupported SNESComposite type");
798:     if (snes->reason < 0) break;

800:     /* Compute the solution update for convergence testing */
801:     VecAXPY(snes->work[0], -1.0, X);
802:     VecScale(snes->work[0], -1.0);

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

807:       if (snes->xl && snes->xu) {
808:         VecNormBegin(X, NORM_2, &xnorm);
809:         VecNormBegin(snes->work[0], NORM_2, &snorm);
810:         SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm);
811:         VecNormEnd(X, NORM_2, &xnorm);
812:         VecNormEnd(snes->work[0], NORM_2, &snorm);
813:       } else {
814:         VecNormBegin(F, NORM_2, &fnorm);
815:         VecNormBegin(X, NORM_2, &xnorm);
816:         VecNormBegin(snes->work[0], NORM_2, &snorm);

818:         VecNormEnd(F, NORM_2, &fnorm);
819:         VecNormEnd(X, NORM_2, &xnorm);
820:         VecNormEnd(snes->work[0], NORM_2, &snorm);
821:       }
822:       SNESCheckFunctionNorm(snes,fnorm);
823:     } else if (normtype == SNES_NORM_ALWAYS) {
824:       VecNormBegin(X, NORM_2, &xnorm);
825:       VecNormBegin(snes->work[0], NORM_2, &snorm);
826:       VecNormEnd(X, NORM_2, &xnorm);
827:       VecNormEnd(snes->work[0], NORM_2, &snorm);
828:     }
829:     /* Monitor convergence */
830:     PetscObjectSAWsTakeAccess((PetscObject)snes);
831:     snes->iter = i+1;
832:     snes->norm = fnorm;
833:     PetscObjectSAWsGrantAccess((PetscObject)snes);
834:     SNESLogConvergenceHistory(snes,snes->norm,0);
835:     SNESMonitor(snes,snes->iter,snes->norm);
836:     /* Test for convergence */
837:     if (normtype == SNES_NORM_ALWAYS) {(*snes->ops->converged)(snes,snes->iter,xnorm,snorm,fnorm,&snes->reason,snes->cnvP);}
838:     if (snes->reason) break;
839:     /* Call general purpose update function */
840:     if (snes->ops->update) {(*snes->ops->update)(snes, snes->iter);}
841:   }
842:   if (normtype == SNES_NORM_ALWAYS) {
843:     if (i == snes->max_its) {
844:       PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);
845:       if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
846:     }
847:   } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS;
848:   return(0);
849: }

851: /* -------------------------------------------------------------------------------------------*/

853: /*MC
854:      SNESCOMPOSITE - Build a preconditioner by composing together several nonlinear solvers

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

860:    Level: intermediate

862:    Concepts: composing solvers

864: .seealso:  SNESCreate(), SNESSetType(), SNESType (for list of available types), SNES,
865:            SNESSHELL, SNESCompositeSetType(), SNESCompositeSpecialSetAlpha(), SNESCompositeAddSNES(),
866:            SNESCompositeGetSNES()

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

872: M*/

874: PETSC_EXTERN PetscErrorCode SNESCreate_Composite(SNES snes)
875: {
877:   SNES_Composite *jac;

880:   PetscNewLog(snes,&jac);

882:   snes->ops->solve           = SNESSolve_Composite;
883:   snes->ops->setup           = SNESSetUp_Composite;
884:   snes->ops->reset           = SNESReset_Composite;
885:   snes->ops->destroy         = SNESDestroy_Composite;
886:   snes->ops->setfromoptions  = SNESSetFromOptions_Composite;
887:   snes->ops->view            = SNESView_Composite;

889:   snes->alwayscomputesfinalresidual = PETSC_FALSE;

891:   snes->data = (void*)jac;
892:   jac->type  = SNES_COMPOSITE_ADDITIVEOPTIMAL;
893:   jac->Fes   = NULL;
894:   jac->Xes   = NULL;
895:   jac->fnorms = NULL;
896:   jac->nsnes = 0;
897:   jac->head  = 0;
898:   jac->stol  = 0.1;
899:   jac->rtol  = 1.1;

901:   jac->h     = NULL;
902:   jac->s     = NULL;
903:   jac->beta  = NULL;
904:   jac->work  = NULL;
905:   jac->rwork = NULL;

907:   PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeSetType_C",SNESCompositeSetType_Composite);
908:   PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeAddSNES_C",SNESCompositeAddSNES_Composite);
909:   PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeGetSNES_C",SNESCompositeGetSNES_Composite);
910:   PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeSetDamping_C",SNESCompositeSetDamping_Composite);
911:   return(0);
912: }