Actual source code: snescomposite.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  2: /*
  3:       Defines a SNES that can consist of a collection of SNESes
  4: */
  5: #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/
  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;

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

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

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

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

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

184: /* approximately solve the overdetermined system:

186:  2*F(x_i)\cdot F(\x_j)\alpha_i = 0
187:  \alpha_i                      = 1

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

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

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

210:   if (snes->normschedule == SNES_NORM_ALWAYS) {
211:     next = jac->head;
212:     SNESSetInitialFunction(next->snes,F);
213:     while (next->next) {
214:       next = next->next;
215:       SNESSetInitialFunction(next->snes,F);
216:     }
217:   }

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

246:   /* all the solutions are collected; combine optimally */
247:   for (i=0;i<jac->n;i++) {
248:     for (j=0;j<i+1;j++) {
249:       VecDotBegin(Fes[i],Fes[j],&jac->h[i + j*jac->n]);
250:     }
251:     VecDotBegin(Fes[i],F,&jac->g[i]);
252:   }

254:   for (i=0;i<jac->n;i++) {
255:     for (j=0;j<i+1;j++) {
256:       VecDotEnd(Fes[i],Fes[j],&jac->h[i + j*jac->n]);
257:       if (i == j) jac->fnorms[i] = PetscSqrtReal(PetscRealPart(jac->h[i + j*jac->n]));
258:     }
259:     VecDotEnd(Fes[i],F,&jac->g[i]);
260:   }

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

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

270:   for (i=0; i<jac->n; i++) {
271:     for (j=0; j<jac->n; j++) {
272:       jac->h[i + j*jac->n] = jac->h[i + j*jac->n] - jac->g[j] - jac->g[i] + ftf;
273:     }
274:     jac->beta[i] = ftf - jac->g[i];
275:   }

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

304:   if (snes->xl && snes->xu) {
305:     SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm);
306:   } else {
307:     VecNorm(F, NORM_2, fnorm);
308:   }

310:   /* take the minimum-normed candidate if it beats the combination by a factor of rtol or the combination has stagnated */
311:   min_fnorm = jac->fnorms[0];
312:   min_i     = 0;
313:   for (i=0; i<jac->n; i++) {
314:     if (jac->fnorms[i] < min_fnorm) {
315:       min_fnorm = jac->fnorms[i];
316:       min_i     = i;
317:     }
318:   }

320:   /* stagnation or divergence restart to the solution of the solver that failed the least */
321:   if (PetscRealPart(total) < jac->stol || min_fnorm*jac->rtol < *fnorm) {
322:     VecCopy(jac->Xes[min_i],X);
323:     VecCopy(jac->Fes[min_i],F);
324:     *fnorm = min_fnorm;
325:   }
326:   return(0);
327: }

331: static PetscErrorCode SNESSetUp_Composite(SNES snes)
332: {
333:   PetscErrorCode     ierr;
334:   DM                 dm;
335:   SNES_Composite     *jac = (SNES_Composite*)snes->data;
336:   SNES_CompositeLink next = jac->head;
337:   PetscInt           n=0,i;
338:   Vec                F;

341:   SNESGetDM(snes,&dm);

343:   if (snes->ops->computevariablebounds) {
344:     /* SNESVI only ever calls computevariablebounds once, so calling it once here is justified */
345:     if (!snes->xl) {VecDuplicate(snes->vec_sol,&snes->xl);}
346:     if (!snes->xu) {VecDuplicate(snes->vec_sol,&snes->xu);}
347:     (*snes->ops->computevariablebounds)(snes,snes->xl,snes->xu);
348:   }

350:   while (next) {
351:     n++;
352:     SNESSetDM(next->snes,dm);
353:     SNESSetApplicationContext(next->snes, snes->user);
354:     if (snes->xl && snes->xu) {
355:       if (snes->ops->computevariablebounds) {
356:         SNESVISetComputeVariableBounds(next->snes, snes->ops->computevariablebounds);
357:       } else {
358:         SNESVISetVariableBounds(next->snes,snes->xl,snes->xu);
359:       }
360:     }

362:     next = next->next;
363:   }
364:   jac->nsnes = n;
365:   SNESGetFunction(snes,&F,NULL,NULL);
366:   if (jac->type == SNES_COMPOSITE_ADDITIVEOPTIMAL) {
367:     VecDuplicateVecs(F,jac->nsnes,&jac->Xes);
368:     PetscMalloc1(n,&jac->Fes);
369:     PetscMalloc1(n,&jac->fnorms);
370:     next = jac->head;
371:     i = 0;
372:     while (next) {
373:       SNESGetFunction(next->snes,&F,NULL,NULL);
374:       jac->Fes[i] = F;
375:       PetscObjectReference((PetscObject)F);
376:       next = next->next;
377:       i++;
378:     }
379:     /* allocate the subspace direct solve area */
380:     jac->nrhs  = 1;
381:     jac->lda   = jac->nsnes;
382:     jac->ldb   = jac->nsnes;
383:     jac->n     = jac->nsnes;

385:     PetscMalloc1(jac->n*jac->n,&jac->h);
386:     PetscMalloc1(jac->n,&jac->beta);
387:     PetscMalloc1(jac->n,&jac->s);
388:     PetscMalloc1(jac->n,&jac->g);
389:     jac->lwork = 12*jac->n;
390: #if PETSC_USE_COMPLEX
391:     PetscMalloc1(jac->lwork,&jac->rwork);
392: #endif
393:     PetscMalloc1(jac->lwork,&jac->work);
394:   }

396:   return(0);
397: }

401: static PetscErrorCode SNESReset_Composite(SNES snes)
402: {
403:   SNES_Composite     *jac = (SNES_Composite*)snes->data;
404:   PetscErrorCode   ierr;
405:   SNES_CompositeLink next = jac->head;

408:   while (next) {
409:     SNESReset(next->snes);
410:     next = next->next;
411:   }
412:   VecDestroy(&jac->Xorig);
413:   if (jac->Xes) {VecDestroyVecs(jac->nsnes,&jac->Xes);}
414:   if (jac->Fes) {VecDestroyVecs(jac->nsnes,&jac->Fes);}
415:   PetscFree(jac->fnorms);
416:   PetscFree(jac->h);
417:   PetscFree(jac->s);
418:   PetscFree(jac->g);
419:   PetscFree(jac->beta);
420:   PetscFree(jac->work);
421:   PetscFree(jac->rwork);
422:   return(0);
423: }

427: static PetscErrorCode SNESDestroy_Composite(SNES snes)
428: {
429:   SNES_Composite     *jac = (SNES_Composite*)snes->data;
430:   PetscErrorCode   ierr;
431:   SNES_CompositeLink next = jac->head,next_tmp;

434:   SNESReset_Composite(snes);
435:   while (next) {
436:     SNESDestroy(&next->snes);
437:     next_tmp = next;
438:     next     = next->next;
439:     PetscFree(next_tmp);
440:   }
441:   PetscFree(snes->data);
442:   return(0);
443: }

447: static PetscErrorCode SNESSetFromOptions_Composite(PetscOptionItems *PetscOptionsObject,SNES snes)
448: {
449:   SNES_Composite     *jac = (SNES_Composite*)snes->data;
450:   PetscErrorCode     ierr;
451:   PetscInt           nmax = 8,i;
452:   SNES_CompositeLink next;
453:   char               *sneses[8];
454:   PetscReal          dmps[8];
455:   PetscBool          flg;

458:   PetscOptionsHead(PetscOptionsObject,"Composite preconditioner options");
459:   PetscOptionsEnum("-snes_composite_type","Type of composition","SNESCompositeSetType",SNESCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&jac->type,&flg);
460:   if (flg) {
461:     SNESCompositeSetType(snes,jac->type);
462:   }
463:   PetscOptionsStringArray("-snes_composite_sneses","List of composite solvers","SNESCompositeAddSNES",sneses,&nmax,&flg);
464:   if (flg) {
465:     for (i=0; i<nmax; i++) {
466:       SNESCompositeAddSNES(snes,sneses[i]);
467:       PetscFree(sneses[i]);   /* deallocate string sneses[i], which is allocated in PetscOptionsStringArray() */
468:     }
469:   }
470:   PetscOptionsRealArray("-snes_composite_damping","Damping of the additive composite solvers","SNESCompositeSetDamping",dmps,&nmax,&flg);
471:   if (flg) {
472:     for (i=0; i<nmax; i++) {
473:       SNESCompositeSetDamping(snes,i,dmps[i]);
474:     }
475:   }
476:   PetscOptionsReal("-snes_composite_stol","Step tolerance for restart on the additive composite solvers","",jac->stol,&jac->stol,NULL);
477:   PetscOptionsReal("-snes_composite_rtol","Residual tolerance for the additive composite solvers","",jac->rtol,&jac->rtol,NULL);
478:   PetscOptionsTail();

480:   next = jac->head;
481:   while (next) {
482:     SNESSetFromOptions(next->snes);
483:     next = next->next;
484:   }
485:   return(0);
486: }

490: static PetscErrorCode SNESView_Composite(SNES snes,PetscViewer viewer)
491: {
492:   SNES_Composite     *jac = (SNES_Composite*)snes->data;
493:   PetscErrorCode   ierr;
494:   SNES_CompositeLink next = jac->head;
495:   PetscBool        iascii;

498:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
499:   if (iascii) {
500:     PetscViewerASCIIPrintf(viewer,"Composite SNES type - %s\n",SNESCompositeTypes[jac->type]);
501:     PetscViewerASCIIPrintf(viewer,"SNESes on composite preconditioner follow\n");
502:     PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
503:   }
504:   if (iascii) {
505:     PetscViewerASCIIPushTab(viewer);
506:   }
507:   while (next) {
508:     SNESView(next->snes,viewer);
509:     next = next->next;
510:   }
511:   if (iascii) {
512:     PetscViewerASCIIPopTab(viewer);
513:     PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
514:   }
515:   return(0);
516: }

518: /* ------------------------------------------------------------------------------*/

522: static PetscErrorCode  SNESCompositeSetType_Composite(SNES snes,SNESCompositeType type)
523: {
524:   SNES_Composite *jac = (SNES_Composite*)snes->data;

527:   jac->type = type;
528:   return(0);
529: }

533: static PetscErrorCode  SNESCompositeAddSNES_Composite(SNES snes,SNESType type)
534: {
535:   SNES_Composite     *jac;
536:   SNES_CompositeLink next,ilink;
537:   PetscErrorCode   ierr;
538:   PetscInt         cnt = 0;
539:   const char       *prefix;
540:   char             newprefix[8];
541:   DM               dm;

544:   PetscNewLog(snes,&ilink);
545:   ilink->next = 0;
546:   SNESCreate(PetscObjectComm((PetscObject)snes),&ilink->snes);
547:   PetscLogObjectParent((PetscObject)snes,(PetscObject)ilink->snes);
548:   SNESGetDM(snes,&dm);
549:   SNESSetDM(ilink->snes,dm);
550:   SNESSetTolerances(ilink->snes,ilink->snes->abstol,ilink->snes->rtol,ilink->snes->stol,1,ilink->snes->max_funcs);
551:   jac  = (SNES_Composite*)snes->data;
552:   next = jac->head;
553:   if (!next) {
554:     jac->head       = ilink;
555:     ilink->previous = NULL;
556:   } else {
557:     cnt++;
558:     while (next->next) {
559:       next = next->next;
560:       cnt++;
561:     }
562:     next->next      = ilink;
563:     ilink->previous = next;
564:   }
565:   SNESGetOptionsPrefix(snes,&prefix);
566:   SNESSetOptionsPrefix(ilink->snes,prefix);
567:   sprintf(newprefix,"sub_%d_",(int)cnt);
568:   SNESAppendOptionsPrefix(ilink->snes,newprefix);
569:   PetscObjectIncrementTabLevel((PetscObject)ilink->snes,(PetscObject)snes,1);
570:   SNESSetType(ilink->snes,type);
571:   SNESSetNormSchedule(ilink->snes, SNES_NORM_FINAL_ONLY);

573:   ilink->dmp = 1.0;
574:   jac->nsnes++;
575:   return(0);
576: }

580: static PetscErrorCode  SNESCompositeGetSNES_Composite(SNES snes,PetscInt n,SNES *subsnes)
581: {
582:   SNES_Composite     *jac;
583:   SNES_CompositeLink next;
584:   PetscInt         i;

587:   jac  = (SNES_Composite*)snes->data;
588:   next = jac->head;
589:   for (i=0; i<n; i++) {
590:     if (!next->next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Not enough SNESes in composite preconditioner");
591:     next = next->next;
592:   }
593:   *subsnes = next->snes;
594:   return(0);
595: }

597: /* -------------------------------------------------------------------------------- */
600: /*@C
601:    SNESCompositeSetType - Sets the type of composite preconditioner.

603:    Logically Collective on SNES

605:    Input Parameter:
606: +  snes - the preconditioner context
607: -  type - SNES_COMPOSITE_ADDITIVE (default), SNES_COMPOSITE_MULTIPLICATIVE

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

612:    Level: Developer

614: .keywords: SNES, set, type, composite preconditioner, additive, multiplicative
615: @*/
616: PetscErrorCode  SNESCompositeSetType(SNES snes,SNESCompositeType type)
617: {

623:   PetscTryMethod(snes,"SNESCompositeSetType_C",(SNES,SNESCompositeType),(snes,type));
624:   return(0);
625: }

629: /*@C
630:    SNESCompositeAddSNES - Adds another SNES to the composite SNES.

632:    Collective on SNES

634:    Input Parameters:
635: +  snes - the preconditioner context
636: -  type - the type of the new preconditioner

638:    Level: Developer

640: .keywords: SNES, composite preconditioner, add
641: @*/
642: PetscErrorCode  SNESCompositeAddSNES(SNES snes,SNESType type)
643: {

648:   PetscTryMethod(snes,"SNESCompositeAddSNES_C",(SNES,SNESType),(snes,type));
649:   return(0);
650: }
653: /*@
654:    SNESCompositeGetSNES - Gets one of the SNES objects in the composite SNES.

656:    Not Collective

658:    Input Parameter:
659: +  snes - the preconditioner context
660: -  n - the number of the snes requested

662:    Output Parameters:
663: .  subsnes - the SNES requested

665:    Level: Developer

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

669: .seealso: SNESCompositeAddSNES()
670: @*/
671: PetscErrorCode  SNESCompositeGetSNES(SNES snes,PetscInt n,SNES *subsnes)
672: {

678:   PetscUseMethod(snes,"SNESCompositeGetSNES_C",(SNES,PetscInt,SNES*),(snes,n,subsnes));
679:   return(0);
680: }

684: /*@
685:    SNESCompositeGetNumber - Get the number of subsolvers in the composite SNES.

687:    Logically Collective on SNES

689:    Input Parameter:
690:    snes - the preconditioner context

692:    Output Parameter:
693:    n - the number of subsolvers

695:    Level: Developer

697: .keywords: SNES, composite preconditioner
698: @*/
699: PetscErrorCode  SNESCompositeGetNumber(SNES snes,PetscInt *n)
700: {
701:   SNES_Composite     *jac;
702:   SNES_CompositeLink next;

705:   jac  = (SNES_Composite*)snes->data;
706:   next = jac->head;

708:   *n = 0;
709:   while (next) {
710:     *n = *n + 1;
711:     next = next->next;
712:   }
713:   return(0);
714: }

718: static PetscErrorCode  SNESCompositeSetDamping_Composite(SNES snes,PetscInt n,PetscReal dmp)
719: {
720:   SNES_Composite     *jac;
721:   SNES_CompositeLink next;
722:   PetscInt         i;

725:   jac  = (SNES_Composite*)snes->data;
726:   next = jac->head;
727:   for (i=0; i<n; i++) {
728:     if (!next->next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Not enough SNESes in composite preconditioner");
729:     next = next->next;
730:   }
731:   next->dmp = dmp;
732:   return(0);
733: }

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

740:    Not Collective

742:    Input Parameter:
743: +  snes - the preconditioner context
744: .  n - the number of the snes requested
745: -  dmp - the damping

747:    Level: Developer

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

751: .seealso: SNESCompositeAddSNES()
752: @*/
753: PetscErrorCode  SNESCompositeSetDamping(SNES snes,PetscInt n,PetscReal dmp)
754: {

759:   PetscUseMethod(snes,"SNESCompositeSetDamping_C",(SNES,PetscInt,PetscReal),(snes,n,dmp));
760:   return(0);
761: }

765: PetscErrorCode SNESSolve_Composite(SNES snes)
766: {
767:   Vec            F;
768:   Vec            X;
769:   Vec            B;
770:   PetscInt       i;
771:   PetscReal      fnorm = 0.0, xnorm = 0.0, snorm = 0.0;
773:   SNESNormSchedule normtype;
774:   SNES_Composite *comp = (SNES_Composite*)snes->data;

777:   X = snes->vec_sol;
778:   F = snes->vec_func;
779:   B = snes->vec_rhs;

781:   PetscObjectSAWsTakeAccess((PetscObject)snes);
782:   snes->iter   = 0;
783:   snes->norm   = 0.;
784:   comp->innerFailures = 0;
785:   PetscObjectSAWsGrantAccess((PetscObject)snes);
786:   SNESSetWorkVecs(snes, 1);
787:   snes->reason = SNES_CONVERGED_ITERATING;
788:   SNESGetNormSchedule(snes, &normtype);
789:   if (normtype == SNES_NORM_ALWAYS || normtype == SNES_NORM_INITIAL_ONLY || normtype == SNES_NORM_INITIAL_FINAL_ONLY) {
790:     if (!snes->vec_func_init_set) {
791:       SNESComputeFunction(snes,X,F);
792:     } else snes->vec_func_init_set = PETSC_FALSE;

794:     if (snes->xl && snes->xu) {
795:       SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm);
796:     } else {
797:       VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
798:     }
799:     SNESCheckFunctionNorm(snes,fnorm);
800:     PetscObjectSAWsTakeAccess((PetscObject)snes);
801:     snes->iter = 0;
802:     snes->norm = fnorm;
803:     PetscObjectSAWsGrantAccess((PetscObject)snes);
804:     SNESLogConvergenceHistory(snes,snes->norm,0);
805:     SNESMonitor(snes,0,snes->norm);

807:     /* test convergence */
808:     (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
809:     if (snes->reason) return(0);
810:   } else {
811:     PetscObjectSAWsGrantAccess((PetscObject)snes);
812:     SNESLogConvergenceHistory(snes,snes->norm,0);
813:     SNESMonitor(snes,0,snes->norm);
814:   }

816:   /* Call general purpose update function */
817:   if (snes->ops->update) {
818:     (*snes->ops->update)(snes, snes->iter);
819:   }

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

826:     if (comp->type == SNES_COMPOSITE_ADDITIVE) {
827:       SNESCompositeApply_Additive(snes,X,B,F,&fnorm);
828:     } else if (comp->type == SNES_COMPOSITE_MULTIPLICATIVE) {
829:       SNESCompositeApply_Multiplicative(snes,X,B,F,&fnorm);
830:     } else if (comp->type == SNES_COMPOSITE_ADDITIVEOPTIMAL) {
831:       SNESCompositeApply_AdditiveOptimal(snes,X,B,F,&fnorm);
832:     } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Unsupported SNESComposite type");
833:     if (snes->reason < 0) break;

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

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

842:       if (snes->xl && snes->xu) {
843:         VecNormBegin(X, NORM_2, &xnorm);
844:         VecNormBegin(snes->work[0], NORM_2, &snorm);
845:         SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm);
846:         VecNormEnd(X, NORM_2, &xnorm);
847:         VecNormEnd(snes->work[0], NORM_2, &snorm);
848:       } else {
849:         VecNormBegin(F, NORM_2, &fnorm);
850:         VecNormBegin(X, NORM_2, &xnorm);
851:         VecNormBegin(snes->work[0], NORM_2, &snorm);

853:         VecNormEnd(F, NORM_2, &fnorm);
854:         VecNormEnd(X, NORM_2, &xnorm);
855:         VecNormEnd(snes->work[0], NORM_2, &snorm);
856:       }
857:       SNESCheckFunctionNorm(snes,fnorm);
858:     } else if (normtype == SNES_NORM_ALWAYS) {
859:       VecNormBegin(X, NORM_2, &xnorm);
860:       VecNormBegin(snes->work[0], NORM_2, &snorm);
861:       VecNormEnd(X, NORM_2, &xnorm);
862:       VecNormEnd(snes->work[0], NORM_2, &snorm);
863:     }
864:     /* Monitor convergence */
865:     PetscObjectSAWsTakeAccess((PetscObject)snes);
866:     snes->iter = i+1;
867:     snes->norm = fnorm;
868:     PetscObjectSAWsGrantAccess((PetscObject)snes);
869:     SNESLogConvergenceHistory(snes,snes->norm,0);
870:     SNESMonitor(snes,snes->iter,snes->norm);
871:     /* Test for convergence */
872:     if (normtype == SNES_NORM_ALWAYS) {(*snes->ops->converged)(snes,snes->iter,xnorm,snorm,fnorm,&snes->reason,snes->cnvP);}
873:     if (snes->reason) break;
874:     /* Call general purpose update function */
875:     if (snes->ops->update) {(*snes->ops->update)(snes, snes->iter);}
876:   }
877:   if (normtype == SNES_NORM_ALWAYS) {
878:     if (i == snes->max_its) {
879:       PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);
880:       if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
881:     }
882:   } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS;
883:   return(0);
884: }

886: /* -------------------------------------------------------------------------------------------*/

888: /*MC
889:      SNESCOMPOSITE - Build a preconditioner by composing together several nonlinear solvers

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

895:    Level: intermediate

897:    Concepts: composing solvers

899: .seealso:  SNESCreate(), SNESSetType(), SNESType (for list of available types), SNES,
900:            SNESSHELL, SNESCompositeSetType(), SNESCompositeSpecialSetAlpha(), SNESCompositeAddSNES(),
901:            SNESCompositeGetSNES()

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

907: M*/

911: PETSC_EXTERN PetscErrorCode SNESCreate_Composite(SNES snes)
912: {
914:   SNES_Composite   *jac;

917:   PetscNewLog(snes,&jac);

919:   snes->ops->solve           = SNESSolve_Composite;
920:   snes->ops->setup           = SNESSetUp_Composite;
921:   snes->ops->reset           = SNESReset_Composite;
922:   snes->ops->destroy         = SNESDestroy_Composite;
923:   snes->ops->setfromoptions  = SNESSetFromOptions_Composite;
924:   snes->ops->view            = SNESView_Composite;

926:   snes->data = (void*)jac;
927:   jac->type  = SNES_COMPOSITE_ADDITIVEOPTIMAL;
928:   jac->Fes   = NULL;
929:   jac->Xes   = NULL;
930:   jac->fnorms = NULL;
931:   jac->nsnes = 0;
932:   jac->head  = 0;
933:   jac->stol  = 0.1;
934:   jac->rtol  = 1.1;

936:   jac->h     = NULL;
937:   jac->s     = NULL;
938:   jac->beta  = NULL;
939:   jac->work  = NULL;
940:   jac->rwork = NULL;

942:   PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeSetType_C",SNESCompositeSetType_Composite);
943:   PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeAddSNES_C",SNESCompositeAddSNES_Composite);
944:   PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeGetSNES_C",SNESCompositeGetSNES_Composite);
945:   PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeSetDamping_C",SNESCompositeSetDamping_Composite);
946:   return(0);
947: }