Actual source code: snes.c

petsc-3.4.4 2014-03-13
  2: #include <petsc-private/snesimpl.h>      /*I "petscsnes.h"  I*/
  3: #include <petscdmshell.h>

  5: PetscBool         SNESRegisterAllCalled = PETSC_FALSE;
  6: PetscFunctionList SNESList              = NULL;

  8: /* Logging support */
  9: PetscClassId  SNES_CLASSID, DMSNES_CLASSID;
 10: PetscLogEvent SNES_Solve, SNES_FunctionEval, SNES_JacobianEval, SNES_GSEval, SNES_NPCSolve;

 14: /*@
 15:    SNESSetErrorIfNotConverged - Causes SNESSolve() to generate an error if the solver has not converged.

 17:    Logically Collective on SNES

 19:    Input Parameters:
 20: +  snes - iterative context obtained from SNESCreate()
 21: -  flg - PETSC_TRUE indicates you want the error generated

 23:    Options database keys:
 24: .  -snes_error_if_not_converged : this takes an optional truth value (0/1/no/yes/true/false)

 26:    Level: intermediate

 28:    Notes:
 29:     Normally PETSc continues if a linear solver fails to converge, you can call SNESGetConvergedReason() after a SNESSolve()
 30:     to determine if it has converged.

 32: .keywords: SNES, set, initial guess, nonzero

 34: .seealso: SNESGetErrorIfNotConverged(), KSPGetErrorIfNotConverged(), KSPSetErrorIFNotConverged()
 35: @*/
 36: PetscErrorCode  SNESSetErrorIfNotConverged(SNES snes,PetscBool flg)
 37: {
 41:   snes->errorifnotconverged = flg;
 42:   return(0);
 43: }

 47: /*@
 48:    SNESGetErrorIfNotConverged - Will SNESSolve() generate an error if the solver does not converge?

 50:    Not Collective

 52:    Input Parameter:
 53: .  snes - iterative context obtained from SNESCreate()

 55:    Output Parameter:
 56: .  flag - PETSC_TRUE if it will generate an error, else PETSC_FALSE

 58:    Level: intermediate

 60: .keywords: SNES, set, initial guess, nonzero

 62: .seealso:  SNESSetErrorIfNotConverged(), KSPGetErrorIfNotConverged(), KSPSetErrorIFNotConverged()
 63: @*/
 64: PetscErrorCode  SNESGetErrorIfNotConverged(SNES snes,PetscBool  *flag)
 65: {
 69:   *flag = snes->errorifnotconverged;
 70:   return(0);
 71: }

 75: /*@
 76:    SNESSetFunctionDomainError - tells SNES that the input vector to your SNESFunction is not
 77:      in the functions domain. For example, negative pressure.

 79:    Logically Collective on SNES

 81:    Input Parameters:
 82: .  snes - the SNES context

 84:    Level: advanced

 86: .keywords: SNES, view

 88: .seealso: SNESCreate(), SNESSetFunction(), SNESFunction
 89: @*/
 90: PetscErrorCode  SNESSetFunctionDomainError(SNES snes)
 91: {
 94:   snes->domainerror = PETSC_TRUE;
 95:   return(0);
 96: }

100: /*@
101:    SNESGetFunctionDomainError - Gets the status of the domain error after a call to SNESComputeFunction;

103:    Logically Collective on SNES

105:    Input Parameters:
106: .  snes - the SNES context

108:    Output Parameters:
109: .  domainerror - Set to PETSC_TRUE if there's a domain error; PETSC_FALSE otherwise.

111:    Level: advanced

113: .keywords: SNES, view

115: .seealso: SNESSetFunctionDomainError(), SNESComputeFunction()
116: @*/
117: PetscErrorCode  SNESGetFunctionDomainError(SNES snes, PetscBool *domainerror)
118: {
122:   *domainerror = snes->domainerror;
123:   return(0);
124: }

128: /*@C
129:   SNESLoad - Loads a SNES that has been stored in binary  with SNESView().

131:   Collective on PetscViewer

133:   Input Parameters:
134: + newdm - the newly loaded SNES, this needs to have been created with SNESCreate() or
135:            some related function before a call to SNESLoad().
136: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen()

138:    Level: intermediate

140:   Notes:
141:    The type is determined by the data in the file, any type set into the SNES before this call is ignored.

143:   Notes for advanced users:
144:   Most users should not need to know the details of the binary storage
145:   format, since SNESLoad() and TSView() completely hide these details.
146:   But for anyone who's interested, the standard binary matrix storage
147:   format is
148: .vb
149:      has not yet been determined
150: .ve

152: .seealso: PetscViewerBinaryOpen(), SNESView(), MatLoad(), VecLoad()
153: @*/
154: PetscErrorCode  SNESLoad(SNES snes, PetscViewer viewer)
155: {
157:   PetscBool      isbinary;
158:   PetscInt       classid;
159:   char           type[256];
160:   KSP            ksp;
161:   DM             dm;
162:   DMSNES         dmsnes;

167:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
168:   if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");

170:   PetscViewerBinaryRead(viewer,&classid,1,PETSC_INT);
171:   if (classid != SNES_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONG,"Not SNES next in file");
172:   PetscViewerBinaryRead(viewer,type,256,PETSC_CHAR);
173:   SNESSetType(snes, type);
174:   if (snes->ops->load) {
175:     (*snes->ops->load)(snes,viewer);
176:   }
177:   SNESGetDM(snes,&dm);
178:   DMGetDMSNES(dm,&dmsnes);
179:   DMSNESLoad(dmsnes,viewer);
180:   SNESGetKSP(snes,&ksp);
181:   KSPLoad(ksp,viewer);
182:   return(0);
183: }

185: #include <petscdraw.h>
186: #if defined(PETSC_HAVE_AMS)
187: #include <petscviewerams.h>
188: #endif
191: /*@C
192:    SNESView - Prints the SNES data structure.

194:    Collective on SNES

196:    Input Parameters:
197: +  SNES - the SNES context
198: -  viewer - visualization context

200:    Options Database Key:
201: .  -snes_view - Calls SNESView() at end of SNESSolve()

203:    Notes:
204:    The available visualization contexts include
205: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
206: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
207:          output where only the first processor opens
208:          the file.  All other processors send their
209:          data to the first processor to print.

211:    The user can open an alternative visualization context with
212:    PetscViewerASCIIOpen() - output to a specified file.

214:    Level: beginner

216: .keywords: SNES, view

218: .seealso: PetscViewerASCIIOpen()
219: @*/
220: PetscErrorCode  SNESView(SNES snes,PetscViewer viewer)
221: {
222:   SNESKSPEW      *kctx;
224:   KSP            ksp;
225:   SNESLineSearch linesearch;
226:   PetscBool      iascii,isstring,isbinary,isdraw;
227:   DMSNES         dmsnes;
228: #if defined(PETSC_HAVE_AMS)
229:   PetscBool      isams;
230: #endif

234:   if (!viewer) {
235:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&viewer);
236:   }

240:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
241:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
242:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
243:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
244: #if defined(PETSC_HAVE_AMS)
245:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERAMS,&isams);
246: #endif
247:   if (iascii) {
248:     PetscObjectPrintClassNamePrefixType((PetscObject)snes,viewer,"SNES Object");
249:     if (!snes->setupcalled) {
250:       PetscViewerASCIIPrintf(viewer,"  SNES has not been set up so information may be incomplete\n");
251:     }
252:     if (snes->ops->view) {
253:       PetscViewerASCIIPushTab(viewer);
254:       (*snes->ops->view)(snes,viewer);
255:       PetscViewerASCIIPopTab(viewer);
256:     }
257:     PetscViewerASCIIPrintf(viewer,"  maximum iterations=%D, maximum function evaluations=%D\n",snes->max_its,snes->max_funcs);
258:     PetscViewerASCIIPrintf(viewer,"  tolerances: relative=%G, absolute=%G, solution=%G\n",snes->rtol,snes->abstol,snes->stol);
259:     PetscViewerASCIIPrintf(viewer,"  total number of linear solver iterations=%D\n",snes->linear_its);
260:     PetscViewerASCIIPrintf(viewer,"  total number of function evaluations=%D\n",snes->nfuncs);
261:     if (snes->gridsequence) {
262:       PetscViewerASCIIPrintf(viewer,"  total number of grid sequence refinements=%D\n",snes->gridsequence);
263:     }
264:     if (snes->ksp_ewconv) {
265:       kctx = (SNESKSPEW*)snes->kspconvctx;
266:       if (kctx) {
267:         PetscViewerASCIIPrintf(viewer,"  Eisenstat-Walker computation of KSP relative tolerance (version %D)\n",kctx->version);
268:         PetscViewerASCIIPrintf(viewer,"    rtol_0=%G, rtol_max=%G, threshold=%G\n",kctx->rtol_0,kctx->rtol_max,kctx->threshold);
269:         PetscViewerASCIIPrintf(viewer,"    gamma=%G, alpha=%G, alpha2=%G\n",kctx->gamma,kctx->alpha,kctx->alpha2);
270:       }
271:     }
272:     if (snes->lagpreconditioner == -1) {
273:       PetscViewerASCIIPrintf(viewer,"  Preconditioned is never rebuilt\n");
274:     } else if (snes->lagpreconditioner > 1) {
275:       PetscViewerASCIIPrintf(viewer,"  Preconditioned is rebuilt every %D new Jacobians\n",snes->lagpreconditioner);
276:     }
277:     if (snes->lagjacobian == -1) {
278:       PetscViewerASCIIPrintf(viewer,"  Jacobian is never rebuilt\n");
279:     } else if (snes->lagjacobian > 1) {
280:       PetscViewerASCIIPrintf(viewer,"  Jacobian is rebuilt every %D SNES iterations\n",snes->lagjacobian);
281:     }
282:   } else if (isstring) {
283:     const char *type;
284:     SNESGetType(snes,&type);
285:     PetscViewerStringSPrintf(viewer," %-3.3s",type);
286:   } else if (isbinary) {
287:     PetscInt    classid = SNES_FILE_CLASSID;
288:     MPI_Comm    comm;
289:     PetscMPIInt rank;
290:     char        type[256];

292:     PetscObjectGetComm((PetscObject)snes,&comm);
293:     MPI_Comm_rank(comm,&rank);
294:     if (!rank) {
295:       PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);
296:       PetscStrncpy(type,((PetscObject)snes)->type_name,256);
297:       PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);
298:     }
299:     if (snes->ops->view) {
300:       (*snes->ops->view)(snes,viewer);
301:     }
302:   } else if (isdraw) {
303:     PetscDraw draw;
304:     char      str[36];
305:     PetscReal x,y,bottom,h;

307:     PetscViewerDrawGetDraw(viewer,0,&draw);
308:     PetscDrawGetCurrentPoint(draw,&x,&y);
309:     PetscStrcpy(str,"SNES: ");
310:     PetscStrcat(str,((PetscObject)snes)->type_name);
311:     PetscDrawBoxedString(draw,x,y,PETSC_DRAW_BLUE,PETSC_DRAW_BLACK,str,NULL,&h);
312:     bottom = y - h;
313:     PetscDrawPushCurrentPoint(draw,x,bottom);
314:     if (snes->ops->view) {
315:       (*snes->ops->view)(snes,viewer);
316:     }
317: #if defined(PETSC_HAVE_AMS)
318:   } else if (isams) {
319:     if (((PetscObject)snes)->amsmem == -1) {
320:       PetscObjectViewAMS((PetscObject)snes,viewer);
321:       PetscStackCallAMS(AMS_Memory_take_access,(((PetscObject)snes)->amsmem));
322:       PetscStackCallAMS(AMS_Memory_add_field,(((PetscObject)snes)->amsmem,"its",&snes->iter,1,AMS_INT,AMS_READ,AMS_COMMON,AMS_REDUCT_UNDEF));
323:       if (!snes->conv_hist) {
324:         SNESSetConvergenceHistory(snes,NULL,NULL,PETSC_DECIDE,PETSC_FALSE);
325:       }
326:       PetscStackCallAMS(AMS_Memory_add_field,(((PetscObject)snes)->amsmem,"conv_hist",snes->conv_hist,10,AMS_DOUBLE,AMS_READ,AMS_COMMON,AMS_REDUCT_UNDEF));
327:       PetscStackCallAMS(AMS_Memory_grant_access,(((PetscObject)snes)->amsmem));
328:     }
329: #endif
330:   }
331:   if (snes->linesearch) {
332:     PetscViewerASCIIPushTab(viewer);
333:     SNESGetLineSearch(snes, &linesearch);
334:     SNESLineSearchView(linesearch, viewer);
335:     PetscViewerASCIIPopTab(viewer);
336:   }
337:   if (snes->pc && snes->usespc) {
338:     PetscViewerASCIIPushTab(viewer);
339:     SNESView(snes->pc, viewer);
340:     PetscViewerASCIIPopTab(viewer);
341:   }
342:   PetscViewerASCIIPushTab(viewer);
343:   DMGetDMSNES(snes->dm,&dmsnes);
344:   DMSNESView(dmsnes, viewer);
345:   PetscViewerASCIIPopTab(viewer);
346:   if (snes->usesksp) {
347:     SNESGetKSP(snes,&ksp);
348:     PetscViewerASCIIPushTab(viewer);
349:     KSPView(ksp,viewer);
350:     PetscViewerASCIIPopTab(viewer);
351:   }
352:   if (isdraw) {
353:     PetscDraw draw;
354:     PetscViewerDrawGetDraw(viewer,0,&draw);
355:     PetscDrawPopCurrentPoint(draw);
356:   }
357:   return(0);
358: }

360: /*
361:   We retain a list of functions that also take SNES command
362:   line options. These are called at the end SNESSetFromOptions()
363: */
364: #define MAXSETFROMOPTIONS 5
365: static PetscInt numberofsetfromoptions;
366: static PetscErrorCode (*othersetfromoptions[MAXSETFROMOPTIONS])(SNES);

370: /*@C
371:   SNESAddOptionsChecker - Adds an additional function to check for SNES options.

373:   Not Collective

375:   Input Parameter:
376: . snescheck - function that checks for options

378:   Level: developer

380: .seealso: SNESSetFromOptions()
381: @*/
382: PetscErrorCode  SNESAddOptionsChecker(PetscErrorCode (*snescheck)(SNES))
383: {
385:   if (numberofsetfromoptions >= MAXSETFROMOPTIONS) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Too many options checkers, only %D allowed", MAXSETFROMOPTIONS);
386:   othersetfromoptions[numberofsetfromoptions++] = snescheck;
387:   return(0);
388: }

390: extern PetscErrorCode  SNESDefaultMatrixFreeCreate2(SNES,Vec,Mat*);

394: static PetscErrorCode SNESSetUpMatrixFree_Private(SNES snes, PetscBool hasOperator, PetscInt version)
395: {
396:   Mat            J;
397:   KSP            ksp;
398:   PC             pc;
399:   PetscBool      match;


405:   if (!snes->vec_func && (snes->jacobian || snes->jacobian_pre)) {
406:     Mat A = snes->jacobian, B = snes->jacobian_pre;
407:     MatGetVecs(A ? A : B, NULL,&snes->vec_func);
408:   }

410:   if (version == 1) {
411:     MatCreateSNESMF(snes,&J);
412:     MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);
413:     MatSetFromOptions(J);
414:   } else if (version == 2) {
415:     if (!snes->vec_func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"SNESSetFunction() must be called first");
416: #if !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128)
417:     SNESDefaultMatrixFreeCreate2(snes,snes->vec_func,&J);
418: #else
419:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "matrix-free operator rutines (version 2)");
420: #endif
421:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "matrix-free operator rutines, only version 1 and 2");

423:   PetscInfo1(snes,"Setting default matrix-free operator routines (version %D)\n", version);
424:   if (hasOperator) {

426:     /* This version replaces the user provided Jacobian matrix with a
427:        matrix-free version but still employs the user-provided preconditioner matrix. */
428:     SNESSetJacobian(snes,J,0,0,0);
429:   } else {
430:     /* This version replaces both the user-provided Jacobian and the user-
431:      provided preconditioner Jacobian with the default matrix free version. */
432:     if ((snes->pcside == PC_LEFT) && snes->pc) {
433:       if (!snes->jacobian){SNESSetJacobian(snes,J,0,0,0);}
434:     } else {
435:       SNESSetJacobian(snes,J,J,MatMFFDComputeJacobian,0);
436:     }
437:     /* Force no preconditioner */
438:     SNESGetKSP(snes,&ksp);
439:     KSPGetPC(ksp,&pc);
440:     PetscObjectTypeCompare((PetscObject)pc,PCSHELL,&match);
441:     if (!match) {
442:       PetscInfo(snes,"Setting default matrix-free preconditioner routines\nThat is no preconditioner is being used\n");
443:       PCSetType(pc,PCNONE);
444:     }
445:   }
446:   MatDestroy(&J);
447:   return(0);
448: }

452: static PetscErrorCode DMRestrictHook_SNESVecSol(DM dmfine,Mat Restrict,Vec Rscale,Mat Inject,DM dmcoarse,void *ctx)
453: {
454:   SNES           snes = (SNES)ctx;
456:   Vec            Xfine,Xfine_named = NULL,Xcoarse;

459:   if (PetscLogPrintInfo) {
460:     PetscInt finelevel,coarselevel,fineclevel,coarseclevel;
461:     DMGetRefineLevel(dmfine,&finelevel);
462:     DMGetCoarsenLevel(dmfine,&fineclevel);
463:     DMGetRefineLevel(dmcoarse,&coarselevel);
464:     DMGetCoarsenLevel(dmcoarse,&coarseclevel);
465:     PetscInfo4(dmfine,"Restricting SNES solution vector from level %D-%D to level %D-%D\n",finelevel,fineclevel,coarselevel,coarseclevel);
466:   }
467:   if (dmfine == snes->dm) Xfine = snes->vec_sol;
468:   else {
469:     DMGetNamedGlobalVector(dmfine,"SNESVecSol",&Xfine_named);
470:     Xfine = Xfine_named;
471:   }
472:   DMGetNamedGlobalVector(dmcoarse,"SNESVecSol",&Xcoarse);
473:   MatRestrict(Restrict,Xfine,Xcoarse);
474:   VecPointwiseMult(Xcoarse,Xcoarse,Rscale);
475:   DMRestoreNamedGlobalVector(dmcoarse,"SNESVecSol",&Xcoarse);
476:   if (Xfine_named) {DMRestoreNamedGlobalVector(dmfine,"SNESVecSol",&Xfine_named);}
477:   return(0);
478: }

482: static PetscErrorCode DMCoarsenHook_SNESVecSol(DM dm,DM dmc,void *ctx)
483: {

487:   DMCoarsenHookAdd(dmc,DMCoarsenHook_SNESVecSol,DMRestrictHook_SNESVecSol,ctx);
488:   return(0);
489: }

493: /* This may be called to rediscretize the operator on levels of linear multigrid. The DM shuffle is so the user can
494:  * safely call SNESGetDM() in their residual evaluation routine. */
495: static PetscErrorCode KSPComputeOperators_SNES(KSP ksp,Mat A,Mat B,MatStructure *mstruct,void *ctx)
496: {
497:   SNES           snes = (SNES)ctx;
499:   Mat            Asave = A,Bsave = B;
500:   Vec            X,Xnamed = NULL;
501:   DM             dmsave;
502:   void           *ctxsave;
503:   PetscErrorCode (*jac)(SNES,Vec,Mat*,Mat*,MatStructure*,void*);

506:   dmsave = snes->dm;
507:   KSPGetDM(ksp,&snes->dm);
508:   if (dmsave == snes->dm) X = snes->vec_sol; /* We are on the finest level */
509:   else {                                     /* We are on a coarser level, this vec was initialized using a DM restrict hook */
510:     DMGetNamedGlobalVector(snes->dm,"SNESVecSol",&Xnamed);
511:     X    = Xnamed;
512:     SNESGetJacobian(snes,NULL,NULL,&jac,&ctxsave);
513:     /* If the DM's don't match up, the MatFDColoring context needed for the jacobian won't match up either -- fixit. */
514:     if (jac == SNESComputeJacobianDefaultColor) {
515:       SNESSetJacobian(snes,NULL,NULL,SNESComputeJacobianDefaultColor,0);
516:     }
517:   }
518:   /* put the previous context back */

520:   SNESComputeJacobian(snes,X,&A,&B,mstruct);
521:   if (snes->dm != dmsave && jac == SNESComputeJacobianDefaultColor) {
522:     SNESSetJacobian(snes,NULL,NULL,jac,ctxsave);
523:   }

525:   if (A != Asave || B != Bsave) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"No support for changing matrices at this time");
526:   if (Xnamed) {
527:     DMRestoreNamedGlobalVector(snes->dm,"SNESVecSol",&Xnamed);
528:   }
529:   snes->dm = dmsave;
530:   return(0);
531: }

535: /*@
536:    SNESSetUpMatrices - ensures that matrices are available for SNES, to be called by SNESSetUp_XXX()

538:    Collective

540:    Input Arguments:
541: .  snes - snes to configure

543:    Level: developer

545: .seealso: SNESSetUp()
546: @*/
547: PetscErrorCode SNESSetUpMatrices(SNES snes)
548: {
550:   DM             dm;
551:   DMSNES         sdm;

554:   SNESGetDM(snes,&dm);
555:   DMGetDMSNES(dm,&sdm);
556:   if (!sdm->ops->computejacobian) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"DMSNES not properly configured");
557:   else if (!snes->jacobian && snes->mf) {
558:     Mat  J;
559:     void *functx;
560:     MatCreateSNESMF(snes,&J);
561:     MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);
562:     MatSetFromOptions(J);
563:     SNESGetFunction(snes,NULL,NULL,&functx);
564:     SNESSetJacobian(snes,J,J,0,0);
565:     MatDestroy(&J);
566:   } else if (snes->mf_operator && !snes->jacobian_pre && !snes->jacobian) {
567:     Mat J,B;
568:     MatCreateSNESMF(snes,&J);
569:     MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);
570:     MatSetFromOptions(J);
571:     DMCreateMatrix(snes->dm,MATAIJ,&B);
572:     /* sdm->computejacobian was already set to reach here */
573:     SNESSetJacobian(snes,J,B,NULL,NULL);
574:     MatDestroy(&J);
575:     MatDestroy(&B);
576:   } else if (!snes->jacobian_pre) {
577:     Mat J,B;
578:     J    = snes->jacobian;
579:     DMCreateMatrix(snes->dm,MATAIJ,&B);
580:     SNESSetJacobian(snes,J ? J : B,B,NULL,NULL);
581:     MatDestroy(&B);
582:   }
583:   {
584:     KSP ksp;
585:     SNESGetKSP(snes,&ksp);
586:     KSPSetComputeOperators(ksp,KSPComputeOperators_SNES,snes);
587:     DMCoarsenHookAdd(snes->dm,DMCoarsenHook_SNESVecSol,DMRestrictHook_SNESVecSol,snes);
588:   }
589:   return(0);
590: }

594: /*@
595:    SNESSetFromOptions - Sets various SNES and KSP parameters from user options.

597:    Collective on SNES

599:    Input Parameter:
600: .  snes - the SNES context

602:    Options Database Keys:
603: +  -snes_type <type> - ls, tr, ngmres, ncg, richardson, qn, vi, fas
604: .  -snes_stol - convergence tolerance in terms of the norm
605:                 of the change in the solution between steps
606: .  -snes_atol <abstol> - absolute tolerance of residual norm
607: .  -snes_rtol <rtol> - relative decrease in tolerance norm from initial
608: .  -snes_max_it <max_it> - maximum number of iterations
609: .  -snes_max_funcs <max_funcs> - maximum number of function evaluations
610: .  -snes_max_fail <max_fail> - maximum number of line search failures allowed before stopping, default is none
611: .  -snes_max_linear_solve_fail - number of linear solver failures before SNESSolve() stops
612: .  -snes_lag_preconditioner <lag> - how often preconditioner is rebuilt (use -1 to never rebuild)
613: .  -snes_lag_jacobian <lag> - how often Jacobian is rebuilt (use -1 to never rebuild)
614: .  -snes_trtol <trtol> - trust region tolerance
615: .  -snes_no_convergence_test - skip convergence test in nonlinear
616:                                solver; hence iterations will continue until max_it
617:                                or some other criterion is reached. Saves expense
618:                                of convergence test
619: .  -snes_monitor <optional filename> - prints residual norm at each iteration. if no
620:                                        filename given prints to stdout
621: .  -snes_monitor_solution - plots solution at each iteration
622: .  -snes_monitor_residual - plots residual (not its norm) at each iteration
623: .  -snes_monitor_solution_update - plots update to solution at each iteration
624: .  -snes_monitor_lg_residualnorm - plots residual norm at each iteration
625: .  -snes_monitor_lg_range - plots residual norm at each iteration
626: .  -snes_fd - use finite differences to compute Jacobian; very slow, only for testing
627: .  -snes_fd_color - use finite differences with coloring to compute Jacobian
628: .  -snes_mf_ksp_monitor - if using matrix-free multiply then print h at each KSP iteration
629: -  -snes_converged_reason - print the reason for convergence/divergence after each solve

631:     Options Database for Eisenstat-Walker method:
632: +  -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence
633: .  -snes_ksp_ew_version ver - version of  Eisenstat-Walker method
634: .  -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0
635: .  -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax
636: .  -snes_ksp_ew_gamma <gamma> - Sets gamma
637: .  -snes_ksp_ew_alpha <alpha> - Sets alpha
638: .  -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2
639: -  -snes_ksp_ew_threshold <threshold> - Sets threshold

641:    Notes:
642:    To see all options, run your program with the -help option or consult
643:    the <A href="../../docs/manual.pdf#nameddest=ch_snes">SNES chapter of the users manual</A>.

645:    Level: beginner

647: .keywords: SNES, nonlinear, set, options, database

649: .seealso: SNESSetOptionsPrefix()
650: @*/
651: PetscErrorCode  SNESSetFromOptions(SNES snes)
652: {
653:   PetscBool      flg,pcset;
654:   PetscInt       i,indx,lag,grids;
655:   MatStructure   matflag;
656:   const char     *deft        = SNESNEWTONLS;
657:   const char     *convtests[] = {"default","skip"};
658:   SNESKSPEW      *kctx        = NULL;
659:   char           type[256], monfilename[PETSC_MAX_PATH_LEN];
660:   PetscViewer    monviewer;
662:   PCSide         pcside;
663:   const char     *optionsprefix;

667:   if (!SNESRegisterAllCalled) {SNESRegisterAll();}
668:   PetscObjectOptionsBegin((PetscObject)snes);
669:   if (((PetscObject)snes)->type_name) deft = ((PetscObject)snes)->type_name;
670:   PetscOptionsList("-snes_type","Nonlinear solver method","SNESSetType",SNESList,deft,type,256,&flg);
671:   if (flg) {
672:     SNESSetType(snes,type);
673:   } else if (!((PetscObject)snes)->type_name) {
674:     SNESSetType(snes,deft);
675:   }
676:   /* not used here, but called so will go into help messaage */
677:   PetscOptionsName("-snes_view","Print detailed information on solver used","SNESView",0);

679:   PetscOptionsReal("-snes_stol","Stop if step length less than","SNESSetTolerances",snes->stol,&snes->stol,0);
680:   PetscOptionsReal("-snes_atol","Stop if function norm less than","SNESSetTolerances",snes->abstol,&snes->abstol,0);

682:   PetscOptionsReal("-snes_rtol","Stop if decrease in function norm less than","SNESSetTolerances",snes->rtol,&snes->rtol,0);
683:   PetscOptionsInt("-snes_max_it","Maximum iterations","SNESSetTolerances",snes->max_its,&snes->max_its,NULL);
684:   PetscOptionsInt("-snes_max_funcs","Maximum function evaluations","SNESSetTolerances",snes->max_funcs,&snes->max_funcs,NULL);
685:   PetscOptionsInt("-snes_max_fail","Maximum nonlinear step failures","SNESSetMaxNonlinearStepFailures",snes->maxFailures,&snes->maxFailures,NULL);
686:   PetscOptionsInt("-snes_max_linear_solve_fail","Maximum failures in linear solves allowed","SNESSetMaxLinearSolveFailures",snes->maxLinearSolveFailures,&snes->maxLinearSolveFailures,NULL);
687:   PetscOptionsBool("-snes_error_if_not_converged","Generate error if solver does not converge","SNESSetErrorIfNotConverged",snes->errorifnotconverged,&snes->errorifnotconverged,NULL);

689:   PetscOptionsInt("-snes_lag_preconditioner","How often to rebuild preconditioner","SNESSetLagPreconditioner",snes->lagpreconditioner,&lag,&flg);
690:   if (flg) {
691:     SNESSetLagPreconditioner(snes,lag);
692:   }
693:   PetscOptionsInt("-snes_lag_jacobian","How often to rebuild Jacobian","SNESSetLagJacobian",snes->lagjacobian,&lag,&flg);
694:   if (flg) {
695:     SNESSetLagJacobian(snes,lag);
696:   }
697:   PetscOptionsInt("-snes_grid_sequence","Use grid sequencing to generate initial guess","SNESSetGridSequence",snes->gridsequence,&grids,&flg);
698:   if (flg) {
699:     SNESSetGridSequence(snes,grids);
700:   }

702:   PetscOptionsEList("-snes_convergence_test","Convergence test","SNESSetConvergenceTest",convtests,2,"default",&indx,&flg);
703:   if (flg) {
704:     switch (indx) {
705:     case 0: SNESSetConvergenceTest(snes,SNESConvergedDefault,NULL,NULL); break;
706:     case 1: SNESSetConvergenceTest(snes,SNESSkipConverged,NULL,NULL);    break;
707:     }
708:   }

710:   PetscOptionsBool("-snes_converged_reason","Print reason for converged or diverged","SNESSolve",snes->printreason,&snes->printreason,NULL);

712:   PetscOptionsEList("-snes_norm_type","SNES Norm type","SNESSetNormType",SNESNormTypes,5,"function",&indx,&flg);
713:   if (flg) { SNESSetNormType(snes,(SNESNormType)indx); }

715:   kctx = (SNESKSPEW*)snes->kspconvctx;

717:   PetscOptionsBool("-snes_ksp_ew","Use Eisentat-Walker linear system convergence test","SNESKSPSetUseEW",snes->ksp_ewconv,&snes->ksp_ewconv,NULL);

719:   PetscOptionsInt("-snes_ksp_ew_version","Version 1, 2 or 3","SNESKSPSetParametersEW",kctx->version,&kctx->version,0);
720:   PetscOptionsReal("-snes_ksp_ew_rtol0","0 <= rtol0 < 1","SNESKSPSetParametersEW",kctx->rtol_0,&kctx->rtol_0,0);
721:   PetscOptionsReal("-snes_ksp_ew_rtolmax","0 <= rtolmax < 1","SNESKSPSetParametersEW",kctx->rtol_max,&kctx->rtol_max,0);
722:   PetscOptionsReal("-snes_ksp_ew_gamma","0 <= gamma <= 1","SNESKSPSetParametersEW",kctx->gamma,&kctx->gamma,0);
723:   PetscOptionsReal("-snes_ksp_ew_alpha","1 < alpha <= 2","SNESKSPSetParametersEW",kctx->alpha,&kctx->alpha,0);
724:   PetscOptionsReal("-snes_ksp_ew_alpha2","alpha2","SNESKSPSetParametersEW",kctx->alpha2,&kctx->alpha2,0);
725:   PetscOptionsReal("-snes_ksp_ew_threshold","0 < threshold < 1","SNESKSPSetParametersEW",kctx->threshold,&kctx->threshold,0);

727:   flg  = PETSC_FALSE;
728:   PetscOptionsBool("-snes_check_jacobian","Check each Jacobian with a differenced one","SNESUpdateCheckJacobian",flg,&flg,NULL);
729:   if (flg) {
730:     SNESSetUpdate(snes,SNESUpdateCheckJacobian);
731:   }

733:   flg  = PETSC_FALSE;
734:   PetscOptionsBool("-snes_monitor_cancel","Remove all monitors","SNESMonitorCancel",flg,&flg,NULL);
735:   if (flg) {SNESMonitorCancel(snes);}

737:   PetscOptionsString("-snes_monitor","Monitor norm of function","SNESMonitorSet","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);
738:   if (flg) {
739:     PetscViewerASCIIOpen(PetscObjectComm((PetscObject)snes),monfilename,&monviewer);
740:     SNESMonitorSet(snes,SNESMonitorDefault,monviewer,(PetscErrorCode (*)(void**))PetscViewerDestroy);
741:   }

743:   PetscOptionsString("-snes_monitor_range","Monitor range of elements of function","SNESMonitorSet","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);
744:   if (flg) {
745:     SNESMonitorSet(snes,SNESMonitorRange,0,0);
746:   }

748:   PetscOptionsString("-snes_ratiomonitor","Monitor ratios of norms of function","SNESMonitorSetRatio","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);
749:   if (flg) {
750:     PetscViewerASCIIOpen(PetscObjectComm((PetscObject)snes),monfilename,&monviewer);
751:     SNESMonitorSetRatio(snes,monviewer);
752:   }

754:   PetscOptionsString("-snes_monitor_short","Monitor norm of function (fewer digits)","SNESMonitorSet","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);
755:   if (flg) {
756:     PetscViewerASCIIOpen(PetscObjectComm((PetscObject)snes),monfilename,&monviewer);
757:     SNESMonitorSet(snes,SNESMonitorDefaultShort,monviewer,(PetscErrorCode (*)(void**))PetscViewerDestroy);
758:   }

760:   PetscOptionsString("-snes_monitor_python","Use Python function","SNESMonitorSet",0,monfilename,PETSC_MAX_PATH_LEN,&flg);
761:   if (flg) {PetscPythonMonitorSet((PetscObject)snes,monfilename);}

763:   flg  = PETSC_FALSE;
764:   PetscOptionsBool("-snes_monitor_solution","Plot solution at each iteration","SNESMonitorSolution",flg,&flg,NULL);
765:   if (flg) {SNESMonitorSet(snes,SNESMonitorSolution,0,0);}
766:   flg  = PETSC_FALSE;
767:   PetscOptionsBool("-snes_monitor_solution_update","Plot correction at each iteration","SNESMonitorSolutionUpdate",flg,&flg,NULL);
768:   if (flg) {SNESMonitorSet(snes,SNESMonitorSolutionUpdate,0,0);}
769:   flg  = PETSC_FALSE;
770:   PetscOptionsBool("-snes_monitor_residual","Plot residual at each iteration","SNESMonitorResidual",flg,&flg,NULL);
771:   if (flg) {SNESMonitorSet(snes,SNESMonitorResidual,0,0);}
772:   flg  = PETSC_FALSE;
773:   PetscOptionsBool("-snes_monitor_lg_residualnorm","Plot function norm at each iteration","SNESMonitorLGResidualNorm",flg,&flg,NULL);
774:   if (flg) {
775:     PetscDrawLG ctx;

777:     SNESMonitorLGCreate(0,0,PETSC_DECIDE,PETSC_DECIDE,300,300,&ctx);
778:     SNESMonitorSet(snes,SNESMonitorLGResidualNorm,ctx,(PetscErrorCode (*)(void**))SNESMonitorLGDestroy);
779:   }
780:   flg  = PETSC_FALSE;
781:   PetscOptionsBool("-snes_monitor_lg_range","Plot function range at each iteration","SNESMonitorLGRange",flg,&flg,NULL);
782:   if (flg) {
783:     PetscViewer ctx;

785:     PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes),0,0,PETSC_DECIDE,PETSC_DECIDE,300,300,&ctx);
786:     SNESMonitorSet(snes,SNESMonitorLGRange,ctx,(PetscErrorCode (*)(void**))PetscViewerDestroy);
787:   }

789:   flg  = PETSC_FALSE;
790:   PetscOptionsBool("-snes_monitor_jacupdate_spectrum","Print the change in the spectrum of the Jacobian","SNESMonitorJacUpdateSpectrum",flg,&flg,NULL);
791:   if (flg) {SNESMonitorSet(snes,SNESMonitorJacUpdateSpectrum,0,0);}

793:   flg  = PETSC_FALSE;
794:   PetscOptionsBool("-snes_fd","Use finite differences (slow) to compute Jacobian","SNESComputeJacobianDefault",flg,&flg,NULL);
795:   if (flg) {
796:     void *functx;
797:     SNESGetFunction(snes,NULL,NULL,&functx);
798:     SNESSetJacobian(snes,snes->jacobian,snes->jacobian_pre,SNESComputeJacobianDefault,functx);
799:     PetscInfo(snes,"Setting default finite difference Jacobian matrix\n");
800:   }

802:   flg  = PETSC_FALSE;
803:   PetscOptionsBool("-snes_fd_function","Use finite differences (slow) to compute function from user objective","SNESObjectiveComputeFunctionDefaultFD",flg,&flg,NULL);
804:   if (flg) {
805:     SNESSetFunction(snes,NULL,SNESObjectiveComputeFunctionDefaultFD,NULL);
806:   }

808:   flg  = PETSC_FALSE;
809:   PetscOptionsBool("-snes_fd_color","Use finite differences with coloring to compute Jacobian","SNESComputeJacobianDefaultColor",flg,&flg,NULL);
810:   if (flg) {
811:     DM             dm;
812:     DMSNES         sdm;
813:     SNESGetDM(snes,&dm);
814:     DMGetDMSNES(dm,&sdm);
815:     sdm->jacobianctx = NULL;
816:     SNESSetJacobian(snes,snes->jacobian,snes->jacobian_pre,SNESComputeJacobianDefaultColor,0);
817:     PetscInfo(snes,"Setting default finite difference coloring Jacobian matrix\n");
818:   }

820:   flg  = PETSC_FALSE;
821:   PetscOptionsBool("-snes_mf_operator","Use a Matrix-Free Jacobian with user-provided preconditioner matrix","MatCreateSNESMF",PETSC_FALSE,&snes->mf_operator,&flg);
822:   if (flg && snes->mf_operator) {
823:     snes->mf_operator = PETSC_TRUE;
824:     snes->mf          = PETSC_TRUE;
825:   }
826:   flg  = PETSC_FALSE;
827:   PetscOptionsBool("-snes_mf","Use a Matrix-Free Jacobian with no preconditioner matrix","MatCreateSNESMF",PETSC_FALSE,&snes->mf,&flg);
828:   if (!flg && snes->mf_operator) snes->mf = PETSC_TRUE;
829:   PetscOptionsInt("-snes_mf_version","Matrix-Free routines version 1 or 2","None",snes->mf_version,&snes->mf_version,0);

831:   flg  = PETSC_FALSE;
832:   SNESGetPCSide(snes,&pcside);
833:   PetscOptionsEnum("-snes_npc_side","SNES nonlinear preconditioner side","SNESSetPCSide",PCSides,(PetscEnum)pcside,(PetscEnum*)&pcside,&flg);
834:   if (flg) {SNESSetPCSide(snes,pcside);}

836: #if defined(PETSC_HAVE_AMS)
837:   {
838:   PetscBool set;
839:   flg  = PETSC_FALSE;
840:   PetscOptionsBool("-snes_ams_block","Block for AMS memory snooper at end of SNESSolve","PetscObjectAMSBlock",((PetscObject)snes)->amspublishblock,&flg,&set);
841:   if (set) {
842:     PetscObjectAMSSetBlock((PetscObject)snes,flg);
843:   }
844:   }
845: #endif

847:   for (i = 0; i < numberofsetfromoptions; i++) {
848:     (*othersetfromoptions[i])(snes);
849:   }

851:   if (snes->ops->setfromoptions) {
852:     (*snes->ops->setfromoptions)(snes);
853:   }

855:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
856:   PetscObjectProcessOptionsHandlers((PetscObject)snes);
857:   PetscOptionsEnd();

859:   if (!snes->ksp) {SNESGetKSP(snes,&snes->ksp);}
860:   KSPGetOperators(snes->ksp,NULL,NULL,&matflag);
861:   KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,matflag);
862:   KSPSetFromOptions(snes->ksp);

864:   if (!snes->linesearch) {
865:     SNESGetLineSearch(snes, &snes->linesearch);
866:   }
867:   SNESLineSearchSetFromOptions(snes->linesearch);

869:   /* if someone has set the SNES PC type, create it. */
870:   SNESGetOptionsPrefix(snes, &optionsprefix);
871:   PetscOptionsHasName(optionsprefix, "-npc_snes_type", &pcset);
872:   if (pcset && (!snes->pc)) {
873:     SNESGetPC(snes, &snes->pc);
874:   }
875:   return(0);
876: }

880: /*@C
881:    SNESSetComputeApplicationContext - Sets an optional function to compute a user-defined context for
882:    the nonlinear solvers.

884:    Logically Collective on SNES

886:    Input Parameters:
887: +  snes - the SNES context
888: .  compute - function to compute the context
889: -  destroy - function to destroy the context

891:    Level: intermediate

893:    Notes:
894:    This function is currently not available from Fortran.

896: .keywords: SNES, nonlinear, set, application, context

898: .seealso: SNESGetApplicationContext(), SNESSetComputeApplicationContext(), SNESGetApplicationContext()
899: @*/
900: PetscErrorCode  SNESSetComputeApplicationContext(SNES snes,PetscErrorCode (*compute)(SNES,void**),PetscErrorCode (*destroy)(void**))
901: {
904:   snes->ops->usercompute = compute;
905:   snes->ops->userdestroy = destroy;
906:   return(0);
907: }

911: /*@
912:    SNESSetApplicationContext - Sets the optional user-defined context for
913:    the nonlinear solvers.

915:    Logically Collective on SNES

917:    Input Parameters:
918: +  snes - the SNES context
919: -  usrP - optional user context

921:    Level: intermediate

923: .keywords: SNES, nonlinear, set, application, context

925: .seealso: SNESGetApplicationContext()
926: @*/
927: PetscErrorCode  SNESSetApplicationContext(SNES snes,void *usrP)
928: {
930:   KSP            ksp;

934:   SNESGetKSP(snes,&ksp);
935:   KSPSetApplicationContext(ksp,usrP);
936:   snes->user = usrP;
937:   return(0);
938: }

942: /*@
943:    SNESGetApplicationContext - Gets the user-defined context for the
944:    nonlinear solvers.

946:    Not Collective

948:    Input Parameter:
949: .  snes - SNES context

951:    Output Parameter:
952: .  usrP - user context

954:    Level: intermediate

956: .keywords: SNES, nonlinear, get, application, context

958: .seealso: SNESSetApplicationContext()
959: @*/
960: PetscErrorCode  SNESGetApplicationContext(SNES snes,void *usrP)
961: {
964:   *(void**)usrP = snes->user;
965:   return(0);
966: }

970: /*@
971:    SNESGetIterationNumber - Gets the number of nonlinear iterations completed
972:    at this time.

974:    Not Collective

976:    Input Parameter:
977: .  snes - SNES context

979:    Output Parameter:
980: .  iter - iteration number

982:    Notes:
983:    For example, during the computation of iteration 2 this would return 1.

985:    This is useful for using lagged Jacobians (where one does not recompute the
986:    Jacobian at each SNES iteration). For example, the code
987: .vb
988:       SNESGetIterationNumber(snes,&it);
989:       if (!(it % 2)) {
990:         [compute Jacobian here]
991:       }
992: .ve
993:    can be used in your ComputeJacobian() function to cause the Jacobian to be
994:    recomputed every second SNES iteration.

996:    Level: intermediate

998: .keywords: SNES, nonlinear, get, iteration, number,

1000: .seealso:   SNESGetFunctionNorm(), SNESGetLinearSolveIterations()
1001: @*/
1002: PetscErrorCode  SNESGetIterationNumber(SNES snes,PetscInt *iter)
1003: {
1007:   *iter = snes->iter;
1008:   return(0);
1009: }

1013: /*@
1014:    SNESSetIterationNumber - Sets the current iteration number.

1016:    Not Collective

1018:    Input Parameter:
1019: .  snes - SNES context
1020: .  iter - iteration number

1022:    Level: developer

1024: .keywords: SNES, nonlinear, set, iteration, number,

1026: .seealso:   SNESGetFunctionNorm(), SNESGetLinearSolveIterations()
1027: @*/
1028: PetscErrorCode  SNESSetIterationNumber(SNES snes,PetscInt iter)
1029: {

1034:   PetscObjectAMSTakeAccess((PetscObject)snes);
1035:   snes->iter = iter;
1036:   PetscObjectAMSGrantAccess((PetscObject)snes);
1037:   return(0);
1038: }

1042: /*@
1043:    SNESGetFunctionNorm - Gets the norm of the current function that was set
1044:    with SNESSSetFunction().

1046:    Collective on SNES

1048:    Input Parameter:
1049: .  snes - SNES context

1051:    Output Parameter:
1052: .  fnorm - 2-norm of function

1054:    Level: intermediate

1056: .keywords: SNES, nonlinear, get, function, norm

1058: .seealso: SNESGetFunction(), SNESGetIterationNumber(), SNESGetLinearSolveIterations()
1059: @*/
1060: PetscErrorCode  SNESGetFunctionNorm(SNES snes,PetscReal *fnorm)
1061: {
1065:   *fnorm = snes->norm;
1066:   return(0);
1067: }


1072: /*@
1073:    SNESSetFunctionNorm - Sets the 2-norm of the current function computed using VecNorm().

1075:    Collective on SNES

1077:    Input Parameter:
1078: .  snes - SNES context
1079: .  fnorm - 2-norm of function

1081:    Level: developer

1083: .keywords: SNES, nonlinear, set, function, norm

1085: .seealso: SNESSetFunction(), SNESSetIterationNumber(), VecNorm().
1086: @*/
1087: PetscErrorCode  SNESSetFunctionNorm(SNES snes,PetscReal fnorm)
1088: {


1094:   PetscObjectAMSTakeAccess((PetscObject)snes);
1095:   snes->norm = fnorm;
1096:   PetscObjectAMSGrantAccess((PetscObject)snes);
1097:   return(0);
1098: }

1102: /*@
1103:    SNESGetNonlinearStepFailures - Gets the number of unsuccessful steps
1104:    attempted by the nonlinear solver.

1106:    Not Collective

1108:    Input Parameter:
1109: .  snes - SNES context

1111:    Output Parameter:
1112: .  nfails - number of unsuccessful steps attempted

1114:    Notes:
1115:    This counter is reset to zero for each successive call to SNESSolve().

1117:    Level: intermediate

1119: .keywords: SNES, nonlinear, get, number, unsuccessful, steps

1121: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures(),
1122:           SNESSetMaxNonlinearStepFailures(), SNESGetMaxNonlinearStepFailures()
1123: @*/
1124: PetscErrorCode  SNESGetNonlinearStepFailures(SNES snes,PetscInt *nfails)
1125: {
1129:   *nfails = snes->numFailures;
1130:   return(0);
1131: }

1135: /*@
1136:    SNESSetMaxNonlinearStepFailures - Sets the maximum number of unsuccessful steps
1137:    attempted by the nonlinear solver before it gives up.

1139:    Not Collective

1141:    Input Parameters:
1142: +  snes     - SNES context
1143: -  maxFails - maximum of unsuccessful steps

1145:    Level: intermediate

1147: .keywords: SNES, nonlinear, set, maximum, unsuccessful, steps

1149: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures(),
1150:           SNESGetMaxNonlinearStepFailures(), SNESGetNonlinearStepFailures()
1151: @*/
1152: PetscErrorCode  SNESSetMaxNonlinearStepFailures(SNES snes, PetscInt maxFails)
1153: {
1156:   snes->maxFailures = maxFails;
1157:   return(0);
1158: }

1162: /*@
1163:    SNESGetMaxNonlinearStepFailures - Gets the maximum number of unsuccessful steps
1164:    attempted by the nonlinear solver before it gives up.

1166:    Not Collective

1168:    Input Parameter:
1169: .  snes     - SNES context

1171:    Output Parameter:
1172: .  maxFails - maximum of unsuccessful steps

1174:    Level: intermediate

1176: .keywords: SNES, nonlinear, get, maximum, unsuccessful, steps

1178: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures(),
1179:           SNESSetMaxNonlinearStepFailures(), SNESGetNonlinearStepFailures()

1181: @*/
1182: PetscErrorCode  SNESGetMaxNonlinearStepFailures(SNES snes, PetscInt *maxFails)
1183: {
1187:   *maxFails = snes->maxFailures;
1188:   return(0);
1189: }

1193: /*@
1194:    SNESGetNumberFunctionEvals - Gets the number of user provided function evaluations
1195:      done by SNES.

1197:    Not Collective

1199:    Input Parameter:
1200: .  snes     - SNES context

1202:    Output Parameter:
1203: .  nfuncs - number of evaluations

1205:    Level: intermediate

1207: .keywords: SNES, nonlinear, get, maximum, unsuccessful, steps

1209: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures()
1210: @*/
1211: PetscErrorCode  SNESGetNumberFunctionEvals(SNES snes, PetscInt *nfuncs)
1212: {
1216:   *nfuncs = snes->nfuncs;
1217:   return(0);
1218: }

1222: /*@
1223:    SNESGetLinearSolveFailures - Gets the number of failed (non-converged)
1224:    linear solvers.

1226:    Not Collective

1228:    Input Parameter:
1229: .  snes - SNES context

1231:    Output Parameter:
1232: .  nfails - number of failed solves

1234:    Notes:
1235:    This counter is reset to zero for each successive call to SNESSolve().

1237:    Level: intermediate

1239: .keywords: SNES, nonlinear, get, number, unsuccessful, steps

1241: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures()
1242: @*/
1243: PetscErrorCode  SNESGetLinearSolveFailures(SNES snes,PetscInt *nfails)
1244: {
1248:   *nfails = snes->numLinearSolveFailures;
1249:   return(0);
1250: }

1254: /*@
1255:    SNESSetMaxLinearSolveFailures - the number of failed linear solve attempts
1256:    allowed before SNES returns with a diverged reason of SNES_DIVERGED_LINEAR_SOLVE

1258:    Logically Collective on SNES

1260:    Input Parameters:
1261: +  snes     - SNES context
1262: -  maxFails - maximum allowed linear solve failures

1264:    Level: intermediate

1266:    Notes: By default this is 0; that is SNES returns on the first failed linear solve

1268: .keywords: SNES, nonlinear, set, maximum, unsuccessful, steps

1270: .seealso: SNESGetLinearSolveFailures(), SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations()
1271: @*/
1272: PetscErrorCode  SNESSetMaxLinearSolveFailures(SNES snes, PetscInt maxFails)
1273: {
1277:   snes->maxLinearSolveFailures = maxFails;
1278:   return(0);
1279: }

1283: /*@
1284:    SNESGetMaxLinearSolveFailures - gets the maximum number of linear solve failures that
1285:      are allowed before SNES terminates

1287:    Not Collective

1289:    Input Parameter:
1290: .  snes     - SNES context

1292:    Output Parameter:
1293: .  maxFails - maximum of unsuccessful solves allowed

1295:    Level: intermediate

1297:    Notes: By default this is 1; that is SNES returns on the first failed linear solve

1299: .keywords: SNES, nonlinear, get, maximum, unsuccessful, steps

1301: .seealso: SNESGetLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(),
1302: @*/
1303: PetscErrorCode  SNESGetMaxLinearSolveFailures(SNES snes, PetscInt *maxFails)
1304: {
1308:   *maxFails = snes->maxLinearSolveFailures;
1309:   return(0);
1310: }

1314: /*@
1315:    SNESGetLinearSolveIterations - Gets the total number of linear iterations
1316:    used by the nonlinear solver.

1318:    Not Collective

1320:    Input Parameter:
1321: .  snes - SNES context

1323:    Output Parameter:
1324: .  lits - number of linear iterations

1326:    Notes:
1327:    This counter is reset to zero for each successive call to SNESSolve().

1329:    Level: intermediate

1331: .keywords: SNES, nonlinear, get, number, linear, iterations

1333: .seealso:  SNESGetIterationNumber(), SNESGetFunctionNorm(), SNESGetLinearSolveFailures(), SNESGetMaxLinearSolveFailures()
1334: @*/
1335: PetscErrorCode  SNESGetLinearSolveIterations(SNES snes,PetscInt *lits)
1336: {
1340:   *lits = snes->linear_its;
1341:   return(0);
1342: }


1347: /*@
1348:    SNESSetKSP - Sets a KSP context for the SNES object to use

1350:    Not Collective, but the SNES and KSP objects must live on the same MPI_Comm

1352:    Input Parameters:
1353: +  snes - the SNES context
1354: -  ksp - the KSP context

1356:    Notes:
1357:    The SNES object already has its KSP object, you can obtain with SNESGetKSP()
1358:    so this routine is rarely needed.

1360:    The KSP object that is already in the SNES object has its reference count
1361:    decreased by one.

1363:    Level: developer

1365: .keywords: SNES, nonlinear, get, KSP, context

1367: .seealso: KSPGetPC(), SNESCreate(), KSPCreate(), SNESSetKSP()
1368: @*/
1369: PetscErrorCode  SNESSetKSP(SNES snes,KSP ksp)
1370: {

1377:   PetscObjectReference((PetscObject)ksp);
1378:   if (snes->ksp) {PetscObjectDereference((PetscObject)snes->ksp);}
1379:   snes->ksp = ksp;
1380:   return(0);
1381: }

1383: /* -----------------------------------------------------------*/
1386: /*@
1387:    SNESCreate - Creates a nonlinear solver context.

1389:    Collective on MPI_Comm

1391:    Input Parameters:
1392: .  comm - MPI communicator

1394:    Output Parameter:
1395: .  outsnes - the new SNES context

1397:    Options Database Keys:
1398: +   -snes_mf - Activates default matrix-free Jacobian-vector products,
1399:                and no preconditioning matrix
1400: .   -snes_mf_operator - Activates default matrix-free Jacobian-vector
1401:                products, and a user-provided preconditioning matrix
1402:                as set by SNESSetJacobian()
1403: -   -snes_fd - Uses (slow!) finite differences to compute Jacobian

1405:    Level: beginner

1407: .keywords: SNES, nonlinear, create, context

1409: .seealso: SNESSolve(), SNESDestroy(), SNES, SNESSetLagPreconditioner()

1411: @*/
1412: PetscErrorCode  SNESCreate(MPI_Comm comm,SNES *outsnes)
1413: {
1415:   SNES           snes;
1416:   SNESKSPEW      *kctx;

1420:   *outsnes = NULL;
1421: #if !defined(PETSC_USE_DYNAMIC_LIBRARIES)
1422:   SNESInitializePackage();
1423: #endif

1425:   PetscHeaderCreate(snes,_p_SNES,struct _SNESOps,SNES_CLASSID,"SNES","Nonlinear solver","SNES",comm,SNESDestroy,SNESView);

1427:   snes->ops->converged    = SNESConvergedDefault;
1428:   snes->usesksp           = PETSC_TRUE;
1429:   snes->tolerancesset     = PETSC_FALSE;
1430:   snes->max_its           = 50;
1431:   snes->max_funcs         = 10000;
1432:   snes->norm              = 0.0;
1433:   snes->normtype          = SNES_NORM_FUNCTION;
1434:   snes->rtol              = 1.e-8;
1435:   snes->ttol              = 0.0;
1436:   snes->abstol            = 1.e-50;
1437:   snes->stol              = 1.e-8;
1438:   snes->deltatol          = 1.e-12;
1439:   snes->nfuncs            = 0;
1440:   snes->numFailures       = 0;
1441:   snes->maxFailures       = 1;
1442:   snes->linear_its        = 0;
1443:   snes->lagjacobian       = 1;
1444:   snes->lagpreconditioner = 1;
1445:   snes->numbermonitors    = 0;
1446:   snes->data              = 0;
1447:   snes->setupcalled       = PETSC_FALSE;
1448:   snes->ksp_ewconv        = PETSC_FALSE;
1449:   snes->nwork             = 0;
1450:   snes->work              = 0;
1451:   snes->nvwork            = 0;
1452:   snes->vwork             = 0;
1453:   snes->conv_hist_len     = 0;
1454:   snes->conv_hist_max     = 0;
1455:   snes->conv_hist         = NULL;
1456:   snes->conv_hist_its     = NULL;
1457:   snes->conv_hist_reset   = PETSC_TRUE;
1458:   snes->vec_func_init_set = PETSC_FALSE;
1459:   snes->norm_init         = 0.;
1460:   snes->norm_init_set     = PETSC_FALSE;
1461:   snes->reason            = SNES_CONVERGED_ITERATING;
1462:   snes->pcside            = PC_RIGHT;

1464:   snes->mf          = PETSC_FALSE;
1465:   snes->mf_operator = PETSC_FALSE;
1466:   snes->mf_version  = 1;

1468:   snes->numLinearSolveFailures = 0;
1469:   snes->maxLinearSolveFailures = 1;

1471:   /* Create context to compute Eisenstat-Walker relative tolerance for KSP */
1472:   PetscNewLog(snes,SNESKSPEW,&kctx);

1474:   snes->kspconvctx  = (void*)kctx;
1475:   kctx->version     = 2;
1476:   kctx->rtol_0      = .3; /* Eisenstat and Walker suggest rtol_0=.5, but
1477:                              this was too large for some test cases */
1478:   kctx->rtol_last   = 0.0;
1479:   kctx->rtol_max    = .9;
1480:   kctx->gamma       = 1.0;
1481:   kctx->alpha       = .5*(1.0 + PetscSqrtReal(5.0));
1482:   kctx->alpha2      = kctx->alpha;
1483:   kctx->threshold   = .1;
1484:   kctx->lresid_last = 0.0;
1485:   kctx->norm_last   = 0.0;

1487:   *outsnes = snes;
1488:   return(0);
1489: }

1491: /*MC
1492:     SNESFunction - functional form used to convey the nonlinear function to be solved by SNES

1494:      Synopsis:
1495:      #include "petscsnes.h"
1496:      SNESFunction(SNES snes,Vec x,Vec f,void *ctx);

1498:      Input Parameters:
1499: +     snes - the SNES context
1500: .     x    - state at which to evaluate residual
1501: -     ctx     - optional user-defined function context, passed in with SNESSetFunction()

1503:      Output Parameter:
1504: .     f  - vector to put residual (function value)

1506:    Level: intermediate

1508: .seealso:   SNESSetFunction(), SNESGetFunction()
1509: M*/

1513: /*@C
1514:    SNESSetFunction - Sets the function evaluation routine and function
1515:    vector for use by the SNES routines in solving systems of nonlinear
1516:    equations.

1518:    Logically Collective on SNES

1520:    Input Parameters:
1521: +  snes - the SNES context
1522: .  r - vector to store function value
1523: .  SNESFunction - function evaluation routine
1524: -  ctx - [optional] user-defined context for private data for the
1525:          function evaluation routine (may be NULL)

1527:    Notes:
1528:    The Newton-like methods typically solve linear systems of the form
1529: $      f'(x) x = -f(x),
1530:    where f'(x) denotes the Jacobian matrix and f(x) is the function.

1532:    Level: beginner

1534: .keywords: SNES, nonlinear, set, function

1536: .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetPicard(), SNESFunction
1537: @*/
1538: PetscErrorCode  SNESSetFunction(SNES snes,Vec r,PetscErrorCode (*SNESFunction)(SNES,Vec,Vec,void*),void *ctx)
1539: {
1541:   DM             dm;

1545:   if (r) {
1548:     PetscObjectReference((PetscObject)r);
1549:     VecDestroy(&snes->vec_func);

1551:     snes->vec_func = r;
1552:   }
1553:   SNESGetDM(snes,&dm);
1554:   DMSNESSetFunction(dm,SNESFunction,ctx);
1555:   return(0);
1556: }


1561: /*@C
1562:    SNESSetInitialFunction - Sets the function vector to be used as the
1563:    function norm at the initialization of the method.  In some
1564:    instances, the user has precomputed the function before calling
1565:    SNESSolve.  This function allows one to avoid a redundant call
1566:    to SNESComputeFunction in that case.

1568:    Logically Collective on SNES

1570:    Input Parameters:
1571: +  snes - the SNES context
1572: -  f - vector to store function value

1574:    Notes:
1575:    This should not be modified during the solution procedure.

1577:    This is used extensively in the SNESFAS hierarchy and in nonlinear preconditioning.

1579:    Level: developer

1581: .keywords: SNES, nonlinear, set, function

1583: .seealso: SNESSetFunction(), SNESComputeFunction(), SNESSetInitialFunctionNorm()
1584: @*/
1585: PetscErrorCode  SNESSetInitialFunction(SNES snes, Vec f)
1586: {
1588:   Vec            vec_func;

1594:   SNESGetFunction(snes,&vec_func,NULL,NULL);
1595:   VecCopy(f, vec_func);

1597:   snes->vec_func_init_set = PETSC_TRUE;
1598:   return(0);
1599: }


1604: /*@C
1605:    SNESSetInitialFunctionNorm - Sets the function norm to be used as the function norm
1606:    at the initialization of the  method.  In some instances, the user has precomputed
1607:    the function and its norm before calling SNESSolve.  This function allows one to
1608:    avoid a redundant call to SNESComputeFunction() and VecNorm() in that case.

1610:    Logically Collective on SNES

1612:    Input Parameters:
1613: +  snes - the SNES context
1614: -  fnorm - the norm of F as set by SNESSetInitialFunction()

1616:    This is used extensively in the SNESFAS hierarchy and in nonlinear preconditioning.

1618:    Level: developer

1620: .keywords: SNES, nonlinear, set, function, norm

1622: .seealso: SNESSetFunction(), SNESComputeFunction(), SNESSetInitialFunction()
1623: @*/
1624: PetscErrorCode  SNESSetInitialFunctionNorm(SNES snes, PetscReal fnorm)
1625: {
1628:   snes->norm_init     = fnorm;
1629:   snes->norm_init_set = PETSC_TRUE;
1630:   return(0);
1631: }

1635: /*@
1636:    SNESSetNormType - Sets the SNESNormType used in covergence and monitoring
1637:    of the SNES method.

1639:    Logically Collective on SNES

1641:    Input Parameters:
1642: +  snes - the SNES context
1643: -  normtype - the type of the norm used

1645:    Notes:
1646:    Only certain SNES methods support certain SNESNormTypes.  Most require evaluation
1647:    of the nonlinear function and the taking of its norm at every iteration to
1648:    even ensure convergence at all.  However, methods such as custom Gauss-Seidel methods
1649:    (SNESGS) and the like do not require the norm of the function to be computed, and therfore
1650:    may either be monitored for convergence or not.  As these are often used as nonlinear
1651:    preconditioners, monitoring the norm of their error is not a useful enterprise within
1652:    their solution.

1654:    Level: developer

1656: .keywords: SNES, nonlinear, set, function, norm, type

1658: .seealso: SNESGetNormType(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormType
1659: @*/
1660: PetscErrorCode  SNESSetNormType(SNES snes, SNESNormType normtype)
1661: {
1664:   snes->normtype = normtype;
1665:   return(0);
1666: }


1671: /*@
1672:    SNESGetNormType - Gets the SNESNormType used in covergence and monitoring
1673:    of the SNES method.

1675:    Logically Collective on SNES

1677:    Input Parameters:
1678: +  snes - the SNES context
1679: -  normtype - the type of the norm used

1681:    Level: advanced

1683: .keywords: SNES, nonlinear, set, function, norm, type

1685: .seealso: SNESSetNormType(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormType
1686: @*/
1687: PetscErrorCode  SNESGetNormType(SNES snes, SNESNormType *normtype)
1688: {
1691:   *normtype = snes->normtype;
1692:   return(0);
1693: }

1695: /*MC
1696:     SNESGSFunction - function used to convey a Gauss-Seidel sweep on the nonlinear function

1698:      Synopsis:
1699:      #include "petscsnes.h"
1700: $    SNESGSFunction(SNES snes,Vec x,Vec b,void *ctx);

1702: +  X   - solution vector
1703: .  B   - RHS vector
1704: -  ctx - optional user-defined Gauss-Seidel context

1706:    Level: intermediate

1708: .seealso:   SNESSetGS(), SNESGetGS()
1709: M*/

1713: /*@C
1714:    SNESSetGS - Sets the user nonlinear Gauss-Seidel routine for
1715:    use with composed nonlinear solvers.

1717:    Input Parameters:
1718: +  snes   - the SNES context
1719: .  SNESGSFunction - function evaluation routine
1720: -  ctx    - [optional] user-defined context for private data for the
1721:             smoother evaluation routine (may be NULL)

1723:    Notes:
1724:    The GS routines are used by the composed nonlinear solver to generate
1725:     a problem appropriate update to the solution, particularly FAS.

1727:    Level: intermediate

1729: .keywords: SNES, nonlinear, set, Gauss-Seidel

1731: .seealso: SNESGetFunction(), SNESComputeGS()
1732: @*/
1733: PetscErrorCode SNESSetGS(SNES snes,PetscErrorCode (*SNESGSFunction)(SNES,Vec,Vec,void*),void *ctx)
1734: {
1736:   DM             dm;

1740:   SNESGetDM(snes,&dm);
1741:   DMSNESSetGS(dm,SNESGSFunction,ctx);
1742:   return(0);
1743: }

1747: PETSC_EXTERN PetscErrorCode SNESPicardComputeFunction(SNES snes,Vec x,Vec f,void *ctx)
1748: {
1750:   DM             dm;
1751:   DMSNES         sdm;

1754:   SNESGetDM(snes,&dm);
1755:   DMGetDMSNES(dm,&sdm);
1756:   /*  A(x)*x - b(x) */
1757:   if (sdm->ops->computepfunction) {
1758:     (*sdm->ops->computepfunction)(snes,x,f,sdm->pctx);
1759:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetPicard() to provide Picard function.");

1761:   if (sdm->ops->computepjacobian) {
1762:     (*sdm->ops->computepjacobian)(snes,x,&snes->jacobian,&snes->jacobian_pre,&snes->matstruct,sdm->pctx);
1763:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetPicard() to provide Picard matrix.");
1764:   VecScale(f,-1.0);
1765:   MatMultAdd(snes->jacobian,x,f,f);
1766:   return(0);
1767: }

1771: PETSC_EXTERN PetscErrorCode SNESPicardComputeJacobian(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx)
1772: {
1774:   /* the jacobian matrix should be pre-filled in SNESPicardComputeFunction */
1775:   *flag = snes->matstruct;
1776:   return(0);
1777: }

1781: /*@C
1782:    SNESSetPicard - Use SNES to solve the semilinear-system A(x) x = b(x) via a Picard type iteration (Picard linearization)

1784:    Logically Collective on SNES

1786:    Input Parameters:
1787: +  snes - the SNES context
1788: .  r - vector to store function value
1789: .  SNESFunction - function evaluation routine
1790: .  Amat - matrix with which A(x) x - b(x) is to be computed
1791: .  Pmat - matrix from which preconditioner is computed (usually the same as Amat)
1792: .  SNESJacobianFunction  - function to compute matrix value
1793: -  ctx - [optional] user-defined context for private data for the
1794:          function evaluation routine (may be NULL)

1796:    Notes:
1797:     We do not recomemend using this routine. It is far better to provide the nonlinear function F() and some approximation to the Jacobian and use
1798:     an approximate Newton solver. This interface is provided to allow porting/testing a previous Picard based code in PETSc before converting it to approximate Newton.

1800:     One can call SNESSetPicard() or SNESSetFunction() (and possibly SNESSetJacobian()) but cannot call both

1802: $     Solves the equation A(x) x = b(x) via the defect correction algorithm A(x^{n}) (x^{n+1} - x^{n}) = b(x^{n}) - A(x^{n})x^{n}
1803: $     Note that when an exact solver is used this corresponds to the "classic" Picard A(x^{n}) x^{n+1} = b(x^{n}) iteration.

1805:      Run with -snes_mf_operator to solve the system with Newton's method using A(x^{n}) to construct the preconditioner.

1807:    We implement the defect correction form of the Picard iteration because it converges much more generally when inexact linear solvers are used then
1808:    the direct Picard iteration A(x^n) x^{n+1} = b(x^n)

1810:    There is some controversity over the definition of a Picard iteration for nonlinear systems but almost everyone agrees that it involves a linear solve and some
1811:    believe it is the iteration  A(x^{n}) x^{n+1} = b(x^{n}) hence we use the name Picard. If anyone has an authoritative  reference that defines the Picard iteration
1812:    different please contact us at petsc-dev@mcs.anl.gov and we'll have an entirely new argument :-).

1814:    Level: intermediate

1816: .keywords: SNES, nonlinear, set, function

1818: .seealso: SNESGetFunction(), SNESSetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESGetPicard(), SNESLineSearchPreCheckPicard()
1819: @*/
1820: PetscErrorCode  SNESSetPicard(SNES snes,Vec r,PetscErrorCode (*SNESFunction)(SNES,Vec,Vec,void*),Mat Amat, Mat Pmat, PetscErrorCode (*SNESJacobianFunction)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void *ctx)
1821: {
1823:   DM             dm;

1827:   SNESGetDM(snes, &dm);
1828:   DMSNESSetPicard(dm,SNESFunction,SNESJacobianFunction,ctx);
1829:   SNESSetFunction(snes,r,SNESPicardComputeFunction,ctx);
1830:   SNESSetJacobian(snes,Amat,Pmat,SNESPicardComputeJacobian,ctx);
1831:   return(0);
1832: }

1836: /*@C
1837:    SNESGetPicard - Returns the context for the Picard iteration

1839:    Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet.

1841:    Input Parameter:
1842: .  snes - the SNES context

1844:    Output Parameter:
1845: +  r - the function (or NULL)
1846: .  SNESFunction - the function (or NULL)
1847: .  Amat - the matrix used to defined the operation A(x) x - b(x) (or NULL)
1848: .  Pmat  - the matrix from which the preconditioner will be constructed (or NULL)
1849: .  SNESJacobianFunction - the function for matrix evaluation (or NULL)
1850: -  ctx - the function context (or NULL)

1852:    Level: advanced

1854: .keywords: SNES, nonlinear, get, function

1856: .seealso: SNESSetPicard(), SNESGetFunction(), SNESGetJacobian(), SNESGetDM(), SNESFunction, SNESJacobianFunction
1857: @*/
1858: PetscErrorCode  SNESGetPicard(SNES snes,Vec *r,PetscErrorCode (**SNESFunction)(SNES,Vec,Vec,void*),Mat *Amat, Mat *Pmat, PetscErrorCode (**SNESJacobianFunction)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void **ctx)
1859: {
1861:   DM             dm;

1865:   SNESGetFunction(snes,r,NULL,NULL);
1866:   SNESGetJacobian(snes,Amat,Pmat,NULL,NULL);
1867:   SNESGetDM(snes,&dm);
1868:   DMSNESGetPicard(dm,SNESFunction,SNESJacobianFunction,ctx);
1869:   return(0);
1870: }

1874: /*@C
1875:    SNESSetComputeInitialGuess - Sets a routine used to compute an initial guess for the problem

1877:    Logically Collective on SNES

1879:    Input Parameters:
1880: +  snes - the SNES context
1881: .  func - function evaluation routine
1882: -  ctx - [optional] user-defined context for private data for the
1883:          function evaluation routine (may be NULL)

1885:    Calling sequence of func:
1886: $    func (SNES snes,Vec x,void *ctx);

1888: .  f - function vector
1889: -  ctx - optional user-defined function context

1891:    Level: intermediate

1893: .keywords: SNES, nonlinear, set, function

1895: .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian()
1896: @*/
1897: PetscErrorCode  SNESSetComputeInitialGuess(SNES snes,PetscErrorCode (*func)(SNES,Vec,void*),void *ctx)
1898: {
1901:   if (func) snes->ops->computeinitialguess = func;
1902:   if (ctx)  snes->initialguessP            = ctx;
1903:   return(0);
1904: }

1906: /* --------------------------------------------------------------- */
1909: /*@C
1910:    SNESGetRhs - Gets the vector for solving F(x) = rhs. If rhs is not set
1911:    it assumes a zero right hand side.

1913:    Logically Collective on SNES

1915:    Input Parameter:
1916: .  snes - the SNES context

1918:    Output Parameter:
1919: .  rhs - the right hand side vector or NULL if the right hand side vector is null

1921:    Level: intermediate

1923: .keywords: SNES, nonlinear, get, function, right hand side

1925: .seealso: SNESGetSolution(), SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction()
1926: @*/
1927: PetscErrorCode  SNESGetRhs(SNES snes,Vec *rhs)
1928: {
1932:   *rhs = snes->vec_rhs;
1933:   return(0);
1934: }

1938: /*@
1939:    SNESComputeFunction - Calls the function that has been set with SNESSetFunction().

1941:    Collective on SNES

1943:    Input Parameters:
1944: +  snes - the SNES context
1945: -  x - input vector

1947:    Output Parameter:
1948: .  y - function vector, as set by SNESSetFunction()

1950:    Notes:
1951:    SNESComputeFunction() is typically used within nonlinear solvers
1952:    implementations, so most users would not generally call this routine
1953:    themselves.

1955:    Level: developer

1957: .keywords: SNES, nonlinear, compute, function

1959: .seealso: SNESSetFunction(), SNESGetFunction()
1960: @*/
1961: PetscErrorCode  SNESComputeFunction(SNES snes,Vec x,Vec y)
1962: {
1964:   DM             dm;
1965:   DMSNES         sdm;

1973:   VecValidValues(x,2,PETSC_TRUE);

1975:   SNESGetDM(snes,&dm);
1976:   DMGetDMSNES(dm,&sdm);
1977:   if (snes->pc && snes->pcside == PC_LEFT) {
1978:     VecCopy(x,y);
1979:     PetscLogEventBegin(SNES_NPCSolve,snes->pc,x,y,0);
1980:     SNESSolve(snes->pc,snes->vec_rhs,y);
1981:     PetscLogEventEnd(SNES_NPCSolve,snes->pc,x,y,0);
1982:     VecAYPX(y,-1.0,x);
1983:     VecValidValues(y,3,PETSC_FALSE);
1984:     return(0);
1985:   } else if (sdm->ops->computefunction) {
1986:     PetscLogEventBegin(SNES_FunctionEval,snes,x,y,0);
1987:     PetscStackPush("SNES user function");
1988:     (*sdm->ops->computefunction)(snes,x,y,sdm->functionctx);
1989:     PetscStackPop;
1990:     PetscLogEventEnd(SNES_FunctionEval,snes,x,y,0);
1991:   } else if (snes->vec_rhs) {
1992:     MatMult(snes->jacobian, x, y);
1993:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetFunction() or SNESSetDM() before SNESComputeFunction(), likely called from SNESSolve().");
1994:   if (snes->vec_rhs) {
1995:     VecAXPY(y,-1.0,snes->vec_rhs);
1996:   }
1997:   snes->nfuncs++;
1998:   VecValidValues(y,3,PETSC_FALSE);
1999:   return(0);
2000: }

2004: /*@
2005:    SNESComputeGS - Calls the Gauss-Seidel function that has been set with  SNESSetGS().

2007:    Collective on SNES

2009:    Input Parameters:
2010: +  snes - the SNES context
2011: .  x - input vector
2012: -  b - rhs vector

2014:    Output Parameter:
2015: .  x - new solution vector

2017:    Notes:
2018:    SNESComputeGS() is typically used within composed nonlinear solver
2019:    implementations, so most users would not generally call this routine
2020:    themselves.

2022:    Level: developer

2024: .keywords: SNES, nonlinear, compute, function

2026: .seealso: SNESSetGS(), SNESComputeFunction()
2027: @*/
2028: PetscErrorCode  SNESComputeGS(SNES snes,Vec b,Vec x)
2029: {
2031:   DM             dm;
2032:   DMSNES         sdm;

2040:   if (b) {VecValidValues(b,2,PETSC_TRUE);}
2041:   PetscLogEventBegin(SNES_GSEval,snes,x,b,0);
2042:   SNESGetDM(snes,&dm);
2043:   DMGetDMSNES(dm,&sdm);
2044:   if (sdm->ops->computegs) {
2045:     PetscStackPush("SNES user GS");
2046:     (*sdm->ops->computegs)(snes,x,b,sdm->gsctx);
2047:     PetscStackPop;
2048:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetGS() before SNESComputeGS(), likely called from SNESSolve().");
2049:   PetscLogEventEnd(SNES_GSEval,snes,x,b,0);
2050:   VecValidValues(x,3,PETSC_FALSE);
2051:   return(0);
2052: }

2056: /*@
2057:    SNESComputeJacobian - Computes the Jacobian matrix that has been set with SNESSetJacobian().

2059:    Collective on SNES and Mat

2061:    Input Parameters:
2062: +  snes - the SNES context
2063: -  x - input vector

2065:    Output Parameters:
2066: +  A - Jacobian matrix
2067: .  B - optional preconditioning matrix
2068: -  flag - flag indicating matrix structure (one of, SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER)

2070:   Options Database Keys:
2071: +    -snes_lag_preconditioner <lag>
2072: .    -snes_lag_jacobian <lag>
2073: .    -snes_compare_explicit - Compare the computed Jacobian to the finite difference Jacobian and output the differences
2074: .    -snes_compare_explicit_draw  - Compare the computed Jacobian to the finite difference Jacobian and draw the result
2075: .    -snes_compare_explicit_contour  - Compare the computed Jacobian to the finite difference Jacobian and draw a contour plot with the result
2076: .    -snes_compare_operator  - Make the comparison options above use the operator instead of the preconditioning matrix
2077: .    -snes_compare_coloring - Compute the finite differece Jacobian using coloring and display norms of difference
2078: .    -snes_compare_coloring_display - Compute the finite differece Jacobian using coloring and display verbose differences
2079: .    -snes_compare_coloring_threshold - Display only those matrix entries that differ by more than a given threshold
2080: .    -snes_compare_coloring_threshold_atol - Absolute tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold
2081: .    -snes_compare_coloring_threshold_rtol - Relative tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold
2082: .    -snes_compare_coloring_draw - Compute the finite differece Jacobian using coloring and draw differences
2083: -    -snes_compare_coloring_draw_contour - Compute the finite differece Jacobian using coloring and show contours of matrices and differences


2086:    Notes:
2087:    Most users should not need to explicitly call this routine, as it
2088:    is used internally within the nonlinear solvers.

2090:    See KSPSetOperators() for important information about setting the
2091:    flag parameter.

2093:    Level: developer

2095: .keywords: SNES, compute, Jacobian, matrix

2097: .seealso:  SNESSetJacobian(), KSPSetOperators(), MatStructure, SNESSetLagPreconditioner(), SNESSetLagJacobian()
2098: @*/
2099: PetscErrorCode  SNESComputeJacobian(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *flg)
2100: {
2102:   PetscBool      flag;
2103:   DM             dm;
2104:   DMSNES         sdm;

2111:   VecValidValues(X,2,PETSC_TRUE);
2112:   SNESGetDM(snes,&dm);
2113:   DMGetDMSNES(dm,&sdm);

2115:   if (!sdm->ops->computejacobian) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_USER,"Must call SNESSetJacobian(), DMSNESSetJacobian(), DMDASNESSetJacobianLocal(), etc");

2117:   /* make sure that MatAssemblyBegin/End() is called on A matrix if it is matrix free */

2119:   if (snes->lagjacobian == -2) {
2120:     snes->lagjacobian = -1;

2122:     PetscInfo(snes,"Recomputing Jacobian/preconditioner because lag is -2 (means compute Jacobian, but then never again) \n");
2123:   } else if (snes->lagjacobian == -1) {
2124:     *flg = SAME_PRECONDITIONER;
2125:     PetscInfo(snes,"Reusing Jacobian/preconditioner because lag is -1\n");
2126:     PetscObjectTypeCompare((PetscObject)*A,MATMFFD,&flag);
2127:     if (flag) {
2128:       MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
2129:       MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
2130:     }
2131:     return(0);
2132:   } else if (snes->lagjacobian > 1 && snes->iter % snes->lagjacobian) {
2133:     *flg = SAME_PRECONDITIONER;
2134:     PetscInfo2(snes,"Reusing Jacobian/preconditioner because lag is %D and SNES iteration is %D\n",snes->lagjacobian,snes->iter);
2135:     PetscObjectTypeCompare((PetscObject)*A,MATMFFD,&flag);
2136:     if (flag) {
2137:       MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
2138:       MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
2139:     }
2140:     return(0);
2141:   }
2142:   if (snes->pc && snes->pcside == PC_LEFT) {
2143:       MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
2144:       MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
2145:       return(0);
2146:   }

2148:   *flg = DIFFERENT_NONZERO_PATTERN;
2149:   PetscLogEventBegin(SNES_JacobianEval,snes,X,*A,*B);

2151:   PetscStackPush("SNES user Jacobian function");
2152:   (*sdm->ops->computejacobian)(snes,X,A,B,flg,sdm->jacobianctx);
2153:   PetscStackPop;

2155:   PetscLogEventEnd(SNES_JacobianEval,snes,X,*A,*B);

2157:   if (snes->lagpreconditioner == -2) {
2158:     PetscInfo(snes,"Rebuilding preconditioner exactly once since lag is -2\n");

2160:     snes->lagpreconditioner = -1;
2161:   } else if (snes->lagpreconditioner == -1) {
2162:     *flg = SAME_PRECONDITIONER;
2163:     PetscInfo(snes,"Reusing preconditioner because lag is -1\n");
2164:   } else if (snes->lagpreconditioner > 1 && snes->iter % snes->lagpreconditioner) {
2165:     *flg = SAME_PRECONDITIONER;
2166:     PetscInfo2(snes,"Reusing preconditioner because lag is %D and SNES iteration is %D\n",snes->lagpreconditioner,snes->iter);
2167:   }

2169:   /* make sure user returned a correct Jacobian and preconditioner */
2172:   {
2173:     PetscBool flag = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_operator = PETSC_FALSE;
2174:     PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit",&flag,NULL);
2175:     PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw",&flag_draw,NULL);
2176:     PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw_contour",&flag_contour,NULL);
2177:     PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_operator",&flag_operator,NULL);
2178:     if (flag || flag_draw || flag_contour) {
2179:       Mat          Bexp_mine = NULL,Bexp,FDexp;
2180:       MatStructure mstruct;
2181:       PetscViewer  vdraw,vstdout;
2182:       PetscBool    flg;
2183:       if (flag_operator) {
2184:         MatComputeExplicitOperator(*A,&Bexp_mine);
2185:         Bexp = Bexp_mine;
2186:       } else {
2187:         /* See if the preconditioning matrix can be viewed and added directly */
2188:         PetscObjectTypeCompareAny((PetscObject)*B,&flg,MATSEQAIJ,MATMPIAIJ,MATSEQDENSE,MATMPIDENSE,MATSEQBAIJ,MATMPIBAIJ,MATSEQSBAIJ,MATMPIBAIJ,"");
2189:         if (flg) Bexp = *B;
2190:         else {
2191:           /* If the "preconditioning" matrix is itself MATSHELL or some other type without direct support */
2192:           MatComputeExplicitOperator(*B,&Bexp_mine);
2193:           Bexp = Bexp_mine;
2194:         }
2195:       }
2196:       MatConvert(Bexp,MATSAME,MAT_INITIAL_MATRIX,&FDexp);
2197:       SNESComputeJacobianDefault(snes,X,&FDexp,&FDexp,&mstruct,NULL);
2198:       PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&vstdout);
2199:       if (flag_draw || flag_contour) {
2200:         PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes),0,"Explicit Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);
2201:         if (flag_contour) {PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);}
2202:       } else vdraw = NULL;
2203:       PetscViewerASCIIPrintf(vstdout,"Explicit %s\n",flag_operator ? "Jacobian" : "preconditioning Jacobian");
2204:       if (flag) {MatView(Bexp,vstdout);}
2205:       if (vdraw) {MatView(Bexp,vdraw);}
2206:       PetscViewerASCIIPrintf(vstdout,"Finite difference Jacobian\n");
2207:       if (flag) {MatView(FDexp,vstdout);}
2208:       if (vdraw) {MatView(FDexp,vdraw);}
2209:       MatAYPX(FDexp,-1.0,Bexp,SAME_NONZERO_PATTERN);
2210:       PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian\n");
2211:       if (flag) {MatView(FDexp,vstdout);}
2212:       if (vdraw) {              /* Always use contour for the difference */
2213:         PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);
2214:         MatView(FDexp,vdraw);
2215:         PetscViewerPopFormat(vdraw);
2216:       }
2217:       if (flag_contour) {PetscViewerPopFormat(vdraw);}
2218:       PetscViewerDestroy(&vdraw);
2219:       MatDestroy(&Bexp_mine);
2220:       MatDestroy(&FDexp);
2221:     }
2222:   }
2223:   {
2224:     PetscBool flag = PETSC_FALSE,flag_display = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_threshold = PETSC_FALSE;
2225:     PetscReal threshold_atol = PETSC_SQRT_MACHINE_EPSILON,threshold_rtol = 10*PETSC_SQRT_MACHINE_EPSILON;
2226:     PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring",&flag,NULL);
2227:     PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_display",&flag_display,NULL);
2228:     PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw",&flag_draw,NULL);
2229:     PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw_contour",&flag_contour,NULL);
2230:     PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold",&flag_threshold,NULL);
2231:     PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_rtol",&threshold_rtol,NULL);
2232:     PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_atol",&threshold_atol,NULL);
2233:     if (flag || flag_display || flag_draw || flag_contour || flag_threshold) {
2234:       Mat            Bfd;
2235:       MatStructure   mstruct;
2236:       PetscViewer    vdraw,vstdout;
2237:       ISColoring     iscoloring;
2238:       MatFDColoring  matfdcoloring;
2239:       PetscErrorCode (*func)(SNES,Vec,Vec,void*);
2240:       void           *funcctx;
2241:       PetscReal      norm1,norm2,normmax;

2243:       MatDuplicate(*B,MAT_DO_NOT_COPY_VALUES,&Bfd);
2244:       MatGetColoring(Bfd,MATCOLORINGSL,&iscoloring);
2245:       MatFDColoringCreate(Bfd,iscoloring,&matfdcoloring);
2246:       ISColoringDestroy(&iscoloring);

2248:       /* This method of getting the function is currently unreliable since it doesn't work for DM local functions. */
2249:       SNESGetFunction(snes,NULL,&func,&funcctx);
2250:       MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))func,funcctx);
2251:       PetscObjectSetOptionsPrefix((PetscObject)matfdcoloring,((PetscObject)snes)->prefix);
2252:       PetscObjectAppendOptionsPrefix((PetscObject)matfdcoloring,"coloring_");
2253:       MatFDColoringSetFromOptions(matfdcoloring);
2254:       MatFDColoringApply(Bfd,matfdcoloring,X,&mstruct,snes);
2255:       MatFDColoringDestroy(&matfdcoloring);

2257:       PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&vstdout);
2258:       if (flag_draw || flag_contour) {
2259:         PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes),0,"Colored Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);
2260:         if (flag_contour) {PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);}
2261:       } else vdraw = NULL;
2262:       PetscViewerASCIIPrintf(vstdout,"Explicit preconditioning Jacobian\n");
2263:       if (flag_display) {MatView(*B,vstdout);}
2264:       if (vdraw) {MatView(*B,vdraw);}
2265:       PetscViewerASCIIPrintf(vstdout,"Colored Finite difference Jacobian\n");
2266:       if (flag_display) {MatView(Bfd,vstdout);}
2267:       if (vdraw) {MatView(Bfd,vdraw);}
2268:       MatAYPX(Bfd,-1.0,*B,SAME_NONZERO_PATTERN);
2269:       MatNorm(Bfd,NORM_1,&norm1);
2270:       MatNorm(Bfd,NORM_FROBENIUS,&norm2);
2271:       MatNorm(Bfd,NORM_MAX,&normmax);
2272:       PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian, norm1=%G normFrob=%G normmax=%G\n",norm1,norm2,normmax);
2273:       if (flag_display) {MatView(Bfd,vstdout);}
2274:       if (vdraw) {              /* Always use contour for the difference */
2275:         PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);
2276:         MatView(Bfd,vdraw);
2277:         PetscViewerPopFormat(vdraw);
2278:       }
2279:       if (flag_contour) {PetscViewerPopFormat(vdraw);}

2281:       if (flag_threshold) {
2282:         PetscInt bs,rstart,rend,i;
2283:         MatGetBlockSize(*B,&bs);
2284:         MatGetOwnershipRange(*B,&rstart,&rend);
2285:         for (i=rstart; i<rend; i++) {
2286:           const PetscScalar *ba,*ca;
2287:           const PetscInt    *bj,*cj;
2288:           PetscInt          bn,cn,j,maxentrycol = -1,maxdiffcol = -1,maxrdiffcol = -1;
2289:           PetscReal         maxentry = 0,maxdiff = 0,maxrdiff = 0;
2290:           MatGetRow(*B,i,&bn,&bj,&ba);
2291:           MatGetRow(Bfd,i,&cn,&cj,&ca);
2292:           if (bn != cn) SETERRQ(((PetscObject)*A)->comm,PETSC_ERR_PLIB,"Unexpected different nonzero pattern in -snes_compare_coloring_threshold");
2293:           for (j=0; j<bn; j++) {
2294:             PetscReal rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j]));
2295:             if (PetscAbsScalar(ba[j]) > PetscAbs(maxentry)) {
2296:               maxentrycol = bj[j];
2297:               maxentry    = PetscRealPart(ba[j]);
2298:             }
2299:             if (PetscAbsScalar(ca[j]) > PetscAbs(maxdiff)) {
2300:               maxdiffcol = bj[j];
2301:               maxdiff    = PetscRealPart(ca[j]);
2302:             }
2303:             if (rdiff > maxrdiff) {
2304:               maxrdiffcol = bj[j];
2305:               maxrdiff    = rdiff;
2306:             }
2307:           }
2308:           if (maxrdiff > 1) {
2309:             PetscViewerASCIIPrintf(vstdout,"row %D (maxentry=%G at %D, maxdiff=%G at %D, maxrdiff=%G at %D):",i,maxentry,maxentrycol,maxdiff,maxdiffcol,maxrdiff,maxrdiffcol);
2310:             for (j=0; j<bn; j++) {
2311:               PetscReal rdiff;
2312:               rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j]));
2313:               if (rdiff > 1) {
2314:                 PetscViewerASCIIPrintf(vstdout," (%D,%G:%G)",bj[j],PetscRealPart(ba[j]),PetscRealPart(ca[j]));
2315:               }
2316:             }
2317:             PetscViewerASCIIPrintf(vstdout,"\n",i,maxentry,maxdiff,maxrdiff);
2318:           }
2319:           MatRestoreRow(*B,i,&bn,&bj,&ba);
2320:           MatRestoreRow(Bfd,i,&cn,&cj,&ca);
2321:         }
2322:       }
2323:       PetscViewerDestroy(&vdraw);
2324:       MatDestroy(&Bfd);
2325:     }
2326:   }
2327:   return(0);
2328: }

2330: /*MC
2331:     SNESJacobianFunction - function used to convey the nonlinear Jacobian of the function to be solved by SNES

2333:      Synopsis:
2334:      #include "petscsnes.h"
2335: $     SNESJacobianFunction(SNES snes,Vec x,Mat *Amat,Mat *Pmat,int *flag,void *ctx);

2337: +  x - input vector
2338: .  Amat - the matrix that defines the (approximate) Jacobian
2339: .  Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.
2340: .  flag - flag indicating information about the preconditioner matrix
2341:    structure (same as flag in KSPSetOperators()), one of SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER
2342: -  ctx - [optional] user-defined Jacobian context

2344:    Level: intermediate

2346: .seealso:   SNESSetFunction(), SNESGetFunction(), SNESSetJacobian(), SNESGetJacobian()
2347: M*/

2351: /*@C
2352:    SNESSetJacobian - Sets the function to compute Jacobian as well as the
2353:    location to store the matrix.

2355:    Logically Collective on SNES and Mat

2357:    Input Parameters:
2358: +  snes - the SNES context
2359: .  Amat - the matrix that defines the (approximate) Jacobian
2360: .  Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.
2361: .  SNESJacobianFunction - Jacobian evaluation routine (if NULL then SNES retains any previously set value)
2362: -  ctx - [optional] user-defined context for private data for the
2363:          Jacobian evaluation routine (may be NULL) (if NULL then SNES retains any previously set value)

2365:    Notes:
2366:    See KSPSetOperators() for important information about setting the flag
2367:    output parameter in the routine func().  Be sure to read this information!

2369:    The routine func() takes Mat * as the matrix arguments rather than Mat.
2370:    This allows the Jacobian evaluation routine to replace A and/or B with a
2371:    completely new new matrix structure (not just different matrix elements)
2372:    when appropriate, for instance, if the nonzero structure is changing
2373:    throughout the global iterations.

2375:    If the Amat matrix and Pmat matrix are different you must call MatAssemblyBegin/End() on
2376:    each matrix.

2378:    If using SNESComputeJacobianDefaultColor() to assemble a Jacobian, the ctx argument
2379:    must be a MatFDColoring.

2381:    Other defect-correction schemes can be used by computing a different matrix in place of the Jacobian.  One common
2382:    example is to use the "Picard linearization" which only differentiates through the highest order parts of each term.

2384:    Level: beginner

2386: .keywords: SNES, nonlinear, set, Jacobian, matrix

2388: .seealso: KSPSetOperators(), SNESSetFunction(), MatMFFDComputeJacobian(), SNESComputeJacobianDefaultColor(), MatStructure, SNESJacobianFunction, SNESSetPicard()
2389: @*/
2390: PetscErrorCode  SNESSetJacobian(SNES snes,Mat Amat,Mat Pmat,PetscErrorCode (*SNESJacobianFunction)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void *ctx)
2391: {
2393:   DM             dm;

2401:   SNESGetDM(snes,&dm);
2402:   DMSNESSetJacobian(dm,SNESJacobianFunction,ctx);
2403:   if (Amat) {
2404:     PetscObjectReference((PetscObject)Amat);
2405:     MatDestroy(&snes->jacobian);

2407:     snes->jacobian = Amat;
2408:   }
2409:   if (Pmat) {
2410:     PetscObjectReference((PetscObject)Pmat);
2411:     MatDestroy(&snes->jacobian_pre);

2413:     snes->jacobian_pre = Pmat;
2414:   }
2415:   return(0);
2416: }

2420: /*@C
2421:    SNESGetJacobian - Returns the Jacobian matrix and optionally the user
2422:    provided context for evaluating the Jacobian.

2424:    Not Collective, but Mat object will be parallel if SNES object is

2426:    Input Parameter:
2427: .  snes - the nonlinear solver context

2429:    Output Parameters:
2430: +  Amat - location to stash (approximate) Jacobian matrix (or NULL)
2431: .  Pmat - location to stash matrix used to compute the preconditioner (or NULL)
2432: .  SNESJacobianFunction - location to put Jacobian function (or NULL)
2433: -  ctx - location to stash Jacobian ctx (or NULL)

2435:    Level: advanced

2437: .seealso: SNESSetJacobian(), SNESComputeJacobian()
2438: @*/
2439: PetscErrorCode SNESGetJacobian(SNES snes,Mat *Amat,Mat *Pmat,PetscErrorCode (**SNESJacobianFunction)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void **ctx)
2440: {
2442:   DM             dm;
2443:   DMSNES         sdm;

2447:   if (Amat) *Amat = snes->jacobian;
2448:   if (Pmat) *Pmat = snes->jacobian_pre;
2449:   SNESGetDM(snes,&dm);
2450:   DMGetDMSNES(dm,&sdm);
2451:   if (SNESJacobianFunction) *SNESJacobianFunction = sdm->ops->computejacobian;
2452:   if (ctx) *ctx = sdm->jacobianctx;
2453:   return(0);
2454: }

2458: /*@
2459:    SNESSetUp - Sets up the internal data structures for the later use
2460:    of a nonlinear solver.

2462:    Collective on SNES

2464:    Input Parameters:
2465: .  snes - the SNES context

2467:    Notes:
2468:    For basic use of the SNES solvers the user need not explicitly call
2469:    SNESSetUp(), since these actions will automatically occur during
2470:    the call to SNESSolve().  However, if one wishes to control this
2471:    phase separately, SNESSetUp() should be called after SNESCreate()
2472:    and optional routines of the form SNESSetXXX(), but before SNESSolve().

2474:    Level: advanced

2476: .keywords: SNES, nonlinear, setup

2478: .seealso: SNESCreate(), SNESSolve(), SNESDestroy()
2479: @*/
2480: PetscErrorCode  SNESSetUp(SNES snes)
2481: {
2483:   DM             dm;
2484:   DMSNES         sdm;
2485:   SNESLineSearch linesearch, pclinesearch;
2486:   void           *lsprectx,*lspostctx;
2487:   PetscErrorCode (*precheck)(SNESLineSearch,Vec,Vec,PetscBool*,void*);
2488:   PetscErrorCode (*postcheck)(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);
2489:   PetscErrorCode (*func)(SNES,Vec,Vec,void*);
2490:   Vec            f,fpc;
2491:   void           *funcctx;
2492:   PetscErrorCode (*jac)(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
2493:   void           *jacctx,*appctx;
2494:   Mat            A,B;

2498:   if (snes->setupcalled) return(0);

2500:   if (!((PetscObject)snes)->type_name) {
2501:     SNESSetType(snes,SNESNEWTONLS);
2502:   }

2504:   SNESGetFunction(snes,&snes->vec_func,NULL,NULL);
2505:   if (snes->vec_func == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be function vector");
2506:   if (snes->vec_rhs  == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be right hand side vector");

2508:   if (!snes->vec_sol_update /* && snes->vec_sol */) {
2509:     VecDuplicate(snes->vec_sol,&snes->vec_sol_update);
2510:     PetscLogObjectParent(snes,snes->vec_sol_update);
2511:   }

2513:   SNESGetDM(snes,&dm);
2514:   DMShellSetGlobalVector(snes->dm,snes->vec_sol);
2515:   DMGetDMSNES(dm,&sdm);
2516:   if (!sdm->ops->computefunction) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Function never provided to SNES object");
2517:   if (!sdm->ops->computejacobian) {
2518:     DMSNESSetJacobian(dm,SNESComputeJacobianDefaultColor,NULL);
2519:   }
2520:   if (!snes->vec_func) {
2521:     DMCreateGlobalVector(dm,&snes->vec_func);
2522:   }

2524:   if (!snes->ksp) {
2525:     SNESGetKSP(snes, &snes->ksp);
2526:   }

2528:   if (!snes->linesearch) {
2529:     SNESGetLineSearch(snes, &snes->linesearch);
2530:   }

2532:   if (snes->pc && (snes->pcside == PC_LEFT)) {
2533:     snes->mf          = PETSC_TRUE;
2534:     snes->mf_operator = PETSC_FALSE;
2535:   }

2537:   if (snes->mf) {
2538:     SNESSetUpMatrixFree_Private(snes, snes->mf_operator, snes->mf_version);
2539:   }

2541:   if (snes->ops->usercompute && !snes->user) {
2542:     (*snes->ops->usercompute)(snes,(void**)&snes->user);
2543:   }

2545:   if (snes->pc) {
2546:     /* copy the DM over */
2547:     SNESGetDM(snes,&dm);
2548:     SNESSetDM(snes->pc,dm);

2550:     SNESGetFunction(snes,&f,&func,&funcctx);
2551:     VecDuplicate(f,&fpc);
2552:     SNESSetFunction(snes->pc,fpc,func,funcctx);
2553:     SNESGetJacobian(snes,&A,&B,&jac,&jacctx);
2554:     SNESSetJacobian(snes->pc,A,B,jac,jacctx);
2555:     SNESGetApplicationContext(snes,&appctx);
2556:     SNESSetApplicationContext(snes->pc,appctx);
2557:     VecDestroy(&fpc);

2559:     /* copy the function pointers over */
2560:     PetscObjectCopyFortranFunctionPointers((PetscObject)snes,(PetscObject)snes->pc);

2562:     /* default to 1 iteration */
2563:     SNESSetTolerances(snes->pc,0.0,0.0,0.0,1,snes->pc->max_funcs);
2564:     if (snes->pcside==PC_RIGHT) {
2565:       SNESSetNormType(snes->pc,SNES_NORM_FINAL_ONLY);
2566:     } else {
2567:       SNESSetNormType(snes->pc,SNES_NORM_NONE);
2568:     }
2569:     SNESSetFromOptions(snes->pc);

2571:     /* copy the line search context over */
2572:     SNESGetLineSearch(snes,&linesearch);
2573:     SNESGetLineSearch(snes->pc,&pclinesearch);
2574:     SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);
2575:     SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);
2576:     SNESLineSearchSetPreCheck(pclinesearch,precheck,lsprectx);
2577:     SNESLineSearchSetPostCheck(pclinesearch,postcheck,lspostctx);
2578:     PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)pclinesearch);
2579:   }

2581:   if (snes->ops->setup) {
2582:     (*snes->ops->setup)(snes);
2583:   }

2585:   snes->setupcalled = PETSC_TRUE;
2586:   return(0);
2587: }

2591: /*@
2592:    SNESReset - Resets a SNES context to the snessetupcalled = 0 state and removes any allocated Vecs and Mats

2594:    Collective on SNES

2596:    Input Parameter:
2597: .  snes - iterative context obtained from SNESCreate()

2599:    Level: intermediate

2601:    Notes: Also calls the application context destroy routine set with SNESSetComputeApplicationContext()

2603: .keywords: SNES, destroy

2605: .seealso: SNESCreate(), SNESSetUp(), SNESSolve()
2606: @*/
2607: PetscErrorCode  SNESReset(SNES snes)
2608: {

2613:   if (snes->ops->userdestroy && snes->user) {
2614:     (*snes->ops->userdestroy)((void**)&snes->user);
2615:     snes->user = NULL;
2616:   }
2617:   if (snes->pc) {
2618:     SNESReset(snes->pc);
2619:   }

2621:   if (snes->ops->reset) {
2622:     (*snes->ops->reset)(snes);
2623:   }
2624:   if (snes->ksp) {
2625:     KSPReset(snes->ksp);
2626:   }

2628:   if (snes->linesearch) {
2629:     SNESLineSearchReset(snes->linesearch);
2630:   }

2632:   VecDestroy(&snes->vec_rhs);
2633:   VecDestroy(&snes->vec_sol);
2634:   VecDestroy(&snes->vec_sol_update);
2635:   VecDestroy(&snes->vec_func);
2636:   MatDestroy(&snes->jacobian);
2637:   MatDestroy(&snes->jacobian_pre);
2638:   VecDestroyVecs(snes->nwork,&snes->work);
2639:   VecDestroyVecs(snes->nvwork,&snes->vwork);

2641:   snes->nwork       = snes->nvwork = 0;
2642:   snes->setupcalled = PETSC_FALSE;
2643:   return(0);
2644: }

2648: /*@
2649:    SNESDestroy - Destroys the nonlinear solver context that was created
2650:    with SNESCreate().

2652:    Collective on SNES

2654:    Input Parameter:
2655: .  snes - the SNES context

2657:    Level: beginner

2659: .keywords: SNES, nonlinear, destroy

2661: .seealso: SNESCreate(), SNESSolve()
2662: @*/
2663: PetscErrorCode  SNESDestroy(SNES *snes)
2664: {

2668:   if (!*snes) return(0);
2670:   if (--((PetscObject)(*snes))->refct > 0) {*snes = 0; return(0);}

2672:   SNESReset((*snes));
2673:   SNESDestroy(&(*snes)->pc);

2675:   /* if memory was published with AMS then destroy it */
2676:   PetscObjectAMSViewOff((PetscObject)*snes);
2677:   if ((*snes)->ops->destroy) {(*((*snes))->ops->destroy)((*snes));}

2679:   DMDestroy(&(*snes)->dm);
2680:   KSPDestroy(&(*snes)->ksp);
2681:   SNESLineSearchDestroy(&(*snes)->linesearch);

2683:   PetscFree((*snes)->kspconvctx);
2684:   if ((*snes)->ops->convergeddestroy) {
2685:     (*(*snes)->ops->convergeddestroy)((*snes)->cnvP);
2686:   }
2687:   if ((*snes)->conv_malloc) {
2688:     PetscFree((*snes)->conv_hist);
2689:     PetscFree((*snes)->conv_hist_its);
2690:   }
2691:   SNESMonitorCancel((*snes));
2692:   PetscHeaderDestroy(snes);
2693:   return(0);
2694: }

2696: /* ----------- Routines to set solver parameters ---------- */

2700: /*@
2701:    SNESSetLagPreconditioner - Determines when the preconditioner is rebuilt in the nonlinear solve.

2703:    Logically Collective on SNES

2705:    Input Parameters:
2706: +  snes - the SNES context
2707: -  lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time
2708:          the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that

2710:    Options Database Keys:
2711: .    -snes_lag_preconditioner <lag>

2713:    Notes:
2714:    The default is 1
2715:    The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1
2716:    If  -1 is used before the very first nonlinear solve the preconditioner is still built because there is no previous preconditioner to use

2718:    Level: intermediate

2720: .keywords: SNES, nonlinear, set, convergence, tolerances

2722: .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian()

2724: @*/
2725: PetscErrorCode  SNESSetLagPreconditioner(SNES snes,PetscInt lag)
2726: {
2729:   if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater");
2730:   if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0");
2732:   snes->lagpreconditioner = lag;
2733:   return(0);
2734: }

2738: /*@
2739:    SNESSetGridSequence - sets the number of steps of grid sequencing that SNES does

2741:    Logically Collective on SNES

2743:    Input Parameters:
2744: +  snes - the SNES context
2745: -  steps - the number of refinements to do, defaults to 0

2747:    Options Database Keys:
2748: .    -snes_grid_sequence <steps>

2750:    Level: intermediate

2752:    Notes:
2753:    Use SNESGetSolution() to extract the fine grid solution after grid sequencing.

2755: .keywords: SNES, nonlinear, set, convergence, tolerances

2757: .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian()

2759: @*/
2760: PetscErrorCode  SNESSetGridSequence(SNES snes,PetscInt steps)
2761: {
2765:   snes->gridsequence = steps;
2766:   return(0);
2767: }

2771: /*@
2772:    SNESGetLagPreconditioner - Indicates how often the preconditioner is rebuilt

2774:    Not Collective

2776:    Input Parameter:
2777: .  snes - the SNES context

2779:    Output Parameter:
2780: .   lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time
2781:          the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that

2783:    Options Database Keys:
2784: .    -snes_lag_preconditioner <lag>

2786:    Notes:
2787:    The default is 1
2788:    The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1

2790:    Level: intermediate

2792: .keywords: SNES, nonlinear, set, convergence, tolerances

2794: .seealso: SNESSetTrustRegionTolerance(), SNESSetLagPreconditioner()

2796: @*/
2797: PetscErrorCode  SNESGetLagPreconditioner(SNES snes,PetscInt *lag)
2798: {
2801:   *lag = snes->lagpreconditioner;
2802:   return(0);
2803: }

2807: /*@
2808:    SNESSetLagJacobian - Determines when the Jacobian is rebuilt in the nonlinear solve. See SNESSetLagPreconditioner() for determining how
2809:      often the preconditioner is rebuilt.

2811:    Logically Collective on SNES

2813:    Input Parameters:
2814: +  snes - the SNES context
2815: -  lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time
2816:          the Jacobian is built etc. -2 means rebuild at next chance but then never again

2818:    Options Database Keys:
2819: .    -snes_lag_jacobian <lag>

2821:    Notes:
2822:    The default is 1
2823:    The Jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1
2824:    If  -1 is used before the very first nonlinear solve the CODE WILL FAIL! because no Jacobian is used, use -2 to indicate you want it recomputed
2825:    at the next Newton step but never again (unless it is reset to another value)

2827:    Level: intermediate

2829: .keywords: SNES, nonlinear, set, convergence, tolerances

2831: .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagPreconditioner(), SNESGetLagJacobian()

2833: @*/
2834: PetscErrorCode  SNESSetLagJacobian(SNES snes,PetscInt lag)
2835: {
2838:   if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater");
2839:   if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0");
2841:   snes->lagjacobian = lag;
2842:   return(0);
2843: }

2847: /*@
2848:    SNESGetLagJacobian - Indicates how often the Jacobian is rebuilt. See SNESGetLagPreconditioner() to determine when the preconditioner is rebuilt

2850:    Not Collective

2852:    Input Parameter:
2853: .  snes - the SNES context

2855:    Output Parameter:
2856: .   lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time
2857:          the Jacobian is built etc.

2859:    Options Database Keys:
2860: .    -snes_lag_jacobian <lag>

2862:    Notes:
2863:    The default is 1
2864:    The jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1

2866:    Level: intermediate

2868: .keywords: SNES, nonlinear, set, convergence, tolerances

2870: .seealso: SNESSetTrustRegionTolerance(), SNESSetLagJacobian(), SNESSetLagPreconditioner(), SNESGetLagPreconditioner()

2872: @*/
2873: PetscErrorCode  SNESGetLagJacobian(SNES snes,PetscInt *lag)
2874: {
2877:   *lag = snes->lagjacobian;
2878:   return(0);
2879: }

2883: /*@
2884:    SNESSetTolerances - Sets various parameters used in convergence tests.

2886:    Logically Collective on SNES

2888:    Input Parameters:
2889: +  snes - the SNES context
2890: .  abstol - absolute convergence tolerance
2891: .  rtol - relative convergence tolerance
2892: .  stol -  convergence tolerance in terms of the norm of the change in the solution between steps,  || delta x || < stol*|| x ||
2893: .  maxit - maximum number of iterations
2894: -  maxf - maximum number of function evaluations

2896:    Options Database Keys:
2897: +    -snes_atol <abstol> - Sets abstol
2898: .    -snes_rtol <rtol> - Sets rtol
2899: .    -snes_stol <stol> - Sets stol
2900: .    -snes_max_it <maxit> - Sets maxit
2901: -    -snes_max_funcs <maxf> - Sets maxf

2903:    Notes:
2904:    The default maximum number of iterations is 50.
2905:    The default maximum number of function evaluations is 1000.

2907:    Level: intermediate

2909: .keywords: SNES, nonlinear, set, convergence, tolerances

2911: .seealso: SNESSetTrustRegionTolerance()
2912: @*/
2913: PetscErrorCode  SNESSetTolerances(SNES snes,PetscReal abstol,PetscReal rtol,PetscReal stol,PetscInt maxit,PetscInt maxf)
2914: {

2923:   if (abstol != PETSC_DEFAULT) {
2924:     if (abstol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %G must be non-negative",abstol);
2925:     snes->abstol = abstol;
2926:   }
2927:   if (rtol != PETSC_DEFAULT) {
2928:     if (rtol < 0.0 || 1.0 <= rtol) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Relative tolerance %G must be non-negative and less than 1.0",rtol);
2929:     snes->rtol = rtol;
2930:   }
2931:   if (stol != PETSC_DEFAULT) {
2932:     if (stol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Step tolerance %G must be non-negative",stol);
2933:     snes->stol = stol;
2934:   }
2935:   if (maxit != PETSC_DEFAULT) {
2936:     if (maxit < 0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxit);
2937:     snes->max_its = maxit;
2938:   }
2939:   if (maxf != PETSC_DEFAULT) {
2940:     if (maxf < 0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of function evaluations %D must be non-negative",maxf);
2941:     snes->max_funcs = maxf;
2942:   }
2943:   snes->tolerancesset = PETSC_TRUE;
2944:   return(0);
2945: }

2949: /*@
2950:    SNESGetTolerances - Gets various parameters used in convergence tests.

2952:    Not Collective

2954:    Input Parameters:
2955: +  snes - the SNES context
2956: .  atol - absolute convergence tolerance
2957: .  rtol - relative convergence tolerance
2958: .  stol -  convergence tolerance in terms of the norm
2959:            of the change in the solution between steps
2960: .  maxit - maximum number of iterations
2961: -  maxf - maximum number of function evaluations

2963:    Notes:
2964:    The user can specify NULL for any parameter that is not needed.

2966:    Level: intermediate

2968: .keywords: SNES, nonlinear, get, convergence, tolerances

2970: .seealso: SNESSetTolerances()
2971: @*/
2972: PetscErrorCode  SNESGetTolerances(SNES snes,PetscReal *atol,PetscReal *rtol,PetscReal *stol,PetscInt *maxit,PetscInt *maxf)
2973: {
2976:   if (atol)  *atol  = snes->abstol;
2977:   if (rtol)  *rtol  = snes->rtol;
2978:   if (stol)  *stol  = snes->stol;
2979:   if (maxit) *maxit = snes->max_its;
2980:   if (maxf)  *maxf  = snes->max_funcs;
2981:   return(0);
2982: }

2986: /*@
2987:    SNESSetTrustRegionTolerance - Sets the trust region parameter tolerance.

2989:    Logically Collective on SNES

2991:    Input Parameters:
2992: +  snes - the SNES context
2993: -  tol - tolerance

2995:    Options Database Key:
2996: .  -snes_trtol <tol> - Sets tol

2998:    Level: intermediate

3000: .keywords: SNES, nonlinear, set, trust region, tolerance

3002: .seealso: SNESSetTolerances()
3003: @*/
3004: PetscErrorCode  SNESSetTrustRegionTolerance(SNES snes,PetscReal tol)
3005: {
3009:   snes->deltatol = tol;
3010:   return(0);
3011: }

3013: /*
3014:    Duplicate the lg monitors for SNES from KSP; for some reason with
3015:    dynamic libraries things don't work under Sun4 if we just use
3016:    macros instead of functions
3017: */
3020: PetscErrorCode  SNESMonitorLGResidualNorm(SNES snes,PetscInt it,PetscReal norm,void *ctx)
3021: {

3026:   KSPMonitorLGResidualNorm((KSP)snes,it,norm,ctx);
3027:   return(0);
3028: }

3032: PetscErrorCode  SNESMonitorLGCreate(const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *draw)
3033: {

3037:   KSPMonitorLGResidualNormCreate(host,label,x,y,m,n,draw);
3038:   return(0);
3039: }

3043: PetscErrorCode  SNESMonitorLGDestroy(PetscDrawLG *draw)
3044: {

3048:   KSPMonitorLGResidualNormDestroy(draw);
3049:   return(0);
3050: }

3052: extern PetscErrorCode  SNESMonitorRange_Private(SNES,PetscInt,PetscReal*);
3055: PetscErrorCode  SNESMonitorLGRange(SNES snes,PetscInt n,PetscReal rnorm,void *monctx)
3056: {
3057:   PetscDrawLG      lg;
3058:   PetscErrorCode   ierr;
3059:   PetscReal        x,y,per;
3060:   PetscViewer      v = (PetscViewer)monctx;
3061:   static PetscReal prev; /* should be in the context */
3062:   PetscDraw        draw;

3065:   PetscViewerDrawGetDrawLG(v,0,&lg);
3066:   if (!n) {PetscDrawLGReset(lg);}
3067:   PetscDrawLGGetDraw(lg,&draw);
3068:   PetscDrawSetTitle(draw,"Residual norm");
3069:   x    = (PetscReal)n;
3070:   if (rnorm > 0.0) y = log10(rnorm);
3071:   else y = -15.0;
3072:   PetscDrawLGAddPoint(lg,&x,&y);
3073:   if (n < 20 || !(n % 5)) {
3074:     PetscDrawLGDraw(lg);
3075:   }

3077:   PetscViewerDrawGetDrawLG(v,1,&lg);
3078:   if (!n) {PetscDrawLGReset(lg);}
3079:   PetscDrawLGGetDraw(lg,&draw);
3080:   PetscDrawSetTitle(draw,"% elemts > .2*max elemt");
3081:    SNESMonitorRange_Private(snes,n,&per);
3082:   x    = (PetscReal)n;
3083:   y    = 100.0*per;
3084:   PetscDrawLGAddPoint(lg,&x,&y);
3085:   if (n < 20 || !(n % 5)) {
3086:     PetscDrawLGDraw(lg);
3087:   }

3089:   PetscViewerDrawGetDrawLG(v,2,&lg);
3090:   if (!n) {prev = rnorm;PetscDrawLGReset(lg);}
3091:   PetscDrawLGGetDraw(lg,&draw);
3092:   PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm");
3093:   x    = (PetscReal)n;
3094:   y    = (prev - rnorm)/prev;
3095:   PetscDrawLGAddPoint(lg,&x,&y);
3096:   if (n < 20 || !(n % 5)) {
3097:     PetscDrawLGDraw(lg);
3098:   }

3100:   PetscViewerDrawGetDrawLG(v,3,&lg);
3101:   if (!n) {PetscDrawLGReset(lg);}
3102:   PetscDrawLGGetDraw(lg,&draw);
3103:   PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm*(% > .2 max)");
3104:   x    = (PetscReal)n;
3105:   y    = (prev - rnorm)/(prev*per);
3106:   if (n > 2) { /*skip initial crazy value */
3107:     PetscDrawLGAddPoint(lg,&x,&y);
3108:   }
3109:   if (n < 20 || !(n % 5)) {
3110:     PetscDrawLGDraw(lg);
3111:   }
3112:   prev = rnorm;
3113:   return(0);
3114: }

3118: /*@
3119:    SNESMonitor - runs the user provided monitor routines, if they exist

3121:    Collective on SNES

3123:    Input Parameters:
3124: +  snes - nonlinear solver context obtained from SNESCreate()
3125: .  iter - iteration number
3126: -  rnorm - relative norm of the residual

3128:    Notes:
3129:    This routine is called by the SNES implementations.
3130:    It does not typically need to be called by the user.

3132:    Level: developer

3134: .seealso: SNESMonitorSet()
3135: @*/
3136: PetscErrorCode  SNESMonitor(SNES snes,PetscInt iter,PetscReal rnorm)
3137: {
3139:   PetscInt       i,n = snes->numbermonitors;

3142:   for (i=0; i<n; i++) {
3143:     (*snes->monitor[i])(snes,iter,rnorm,snes->monitorcontext[i]);
3144:   }
3145:   return(0);
3146: }

3148: /* ------------ Routines to set performance monitoring options ----------- */

3150: /*MC
3151:     SNESMonitorFunction - functional form passed to SNESMonitorSet() to monitor convergence of nonlinear solver

3153:      Synopsis:
3154:      #include "petscsnes.h"
3155: $    PetscErrorCode SNESMonitorFunction(SNES snes,PetscInt its, PetscReal norm,void *mctx)

3157: +    snes - the SNES context
3158: .    its - iteration number
3159: .    norm - 2-norm function value (may be estimated)
3160: -    mctx - [optional] monitoring context

3162:    Level: advanced

3164: .seealso:   SNESMonitorSet(), SNESMonitorGet()
3165: M*/

3169: /*@C
3170:    SNESMonitorSet - Sets an ADDITIONAL function that is to be used at every
3171:    iteration of the nonlinear solver to display the iteration's
3172:    progress.

3174:    Logically Collective on SNES

3176:    Input Parameters:
3177: +  snes - the SNES context
3178: .  SNESMonitorFunction - monitoring routine
3179: .  mctx - [optional] user-defined context for private data for the
3180:           monitor routine (use NULL if no context is desired)
3181: -  monitordestroy - [optional] routine that frees monitor context
3182:           (may be NULL)

3184:    Options Database Keys:
3185: +    -snes_monitor        - sets SNESMonitorDefault()
3186: .    -snes_monitor_lg_residualnorm    - sets line graph monitor,
3187:                             uses SNESMonitorLGCreate()
3188: -    -snes_monitor_cancel - cancels all monitors that have
3189:                             been hardwired into a code by
3190:                             calls to SNESMonitorSet(), but
3191:                             does not cancel those set via
3192:                             the options database.

3194:    Notes:
3195:    Several different monitoring routines may be set by calling
3196:    SNESMonitorSet() multiple times; all will be called in the
3197:    order in which they were set.

3199:    Fortran notes: Only a single monitor function can be set for each SNES object

3201:    Level: intermediate

3203: .keywords: SNES, nonlinear, set, monitor

3205: .seealso: SNESMonitorDefault(), SNESMonitorCancel(), SNESMonitorFunction
3206: @*/
3207: PetscErrorCode  SNESMonitorSet(SNES snes,PetscErrorCode (*SNESMonitorFunction)(SNES,PetscInt,PetscReal,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**))
3208: {
3209:   PetscInt       i;

3214:   if (snes->numbermonitors >= MAXSNESMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set");
3215:   for (i=0; i<snes->numbermonitors;i++) {
3216:     if (SNESMonitorFunction == snes->monitor[i] && monitordestroy == snes->monitordestroy[i] && mctx == snes->monitorcontext[i]) {
3217:       if (monitordestroy) {
3218:         (*monitordestroy)(&mctx);
3219:       }
3220:       return(0);
3221:     }
3222:   }
3223:   snes->monitor[snes->numbermonitors]          = SNESMonitorFunction;
3224:   snes->monitordestroy[snes->numbermonitors]   = monitordestroy;
3225:   snes->monitorcontext[snes->numbermonitors++] = (void*)mctx;
3226:   return(0);
3227: }

3231: /*@
3232:    SNESMonitorCancel - Clears all the monitor functions for a SNES object.

3234:    Logically Collective on SNES

3236:    Input Parameters:
3237: .  snes - the SNES context

3239:    Options Database Key:
3240: .  -snes_monitor_cancel - cancels all monitors that have been hardwired
3241:     into a code by calls to SNESMonitorSet(), but does not cancel those
3242:     set via the options database

3244:    Notes:
3245:    There is no way to clear one specific monitor from a SNES object.

3247:    Level: intermediate

3249: .keywords: SNES, nonlinear, set, monitor

3251: .seealso: SNESMonitorDefault(), SNESMonitorSet()
3252: @*/
3253: PetscErrorCode  SNESMonitorCancel(SNES snes)
3254: {
3256:   PetscInt       i;

3260:   for (i=0; i<snes->numbermonitors; i++) {
3261:     if (snes->monitordestroy[i]) {
3262:       (*snes->monitordestroy[i])(&snes->monitorcontext[i]);
3263:     }
3264:   }
3265:   snes->numbermonitors = 0;
3266:   return(0);
3267: }

3269: /*MC
3270:     SNESConvergenceTestFunction - functional form used for testing of convergence of nonlinear solver

3272:      Synopsis:
3273:      #include "petscsnes.h"
3274: $     PetscErrorCode SNESConvergenceTest(SNES snes,PetscInt it,PetscReal xnorm,PetscReal gnorm,PetscReal f,SNESConvergedReason *reason,void *cctx)

3276: +    snes - the SNES context
3277: .    it - current iteration (0 is the first and is before any Newton step)
3278: .    cctx - [optional] convergence context
3279: .    reason - reason for convergence/divergence
3280: .    xnorm - 2-norm of current iterate
3281: .    gnorm - 2-norm of current step
3282: -    f - 2-norm of function

3284:    Level: intermediate

3286: .seealso:   SNESSetConvergenceTest(), SNESGetConvergenceTest()
3287: M*/

3291: /*@C
3292:    SNESSetConvergenceTest - Sets the function that is to be used
3293:    to test for convergence of the nonlinear iterative solution.

3295:    Logically Collective on SNES

3297:    Input Parameters:
3298: +  snes - the SNES context
3299: .  SNESConvergenceTestFunction - routine to test for convergence
3300: .  cctx - [optional] context for private data for the convergence routine  (may be NULL)
3301: -  destroy - [optional] destructor for the context (may be NULL; NULL_FUNCTION in Fortran)

3303:    Level: advanced

3305: .keywords: SNES, nonlinear, set, convergence, test

3307: .seealso: SNESConvergedDefault(), SNESSkipConverged(), SNESConvergenceTestFunction
3308: @*/
3309: PetscErrorCode  SNESSetConvergenceTest(SNES snes,PetscErrorCode (*SNESConvergenceTestFunction)(SNES,PetscInt,PetscReal,PetscReal,PetscReal,SNESConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*))
3310: {

3315:   if (!SNESConvergenceTestFunction) SNESConvergenceTestFunction = SNESSkipConverged;
3316:   if (snes->ops->convergeddestroy) {
3317:     (*snes->ops->convergeddestroy)(snes->cnvP);
3318:   }
3319:   snes->ops->converged        = SNESConvergenceTestFunction;
3320:   snes->ops->convergeddestroy = destroy;
3321:   snes->cnvP                  = cctx;
3322:   return(0);
3323: }

3327: /*@
3328:    SNESGetConvergedReason - Gets the reason the SNES iteration was stopped.

3330:    Not Collective

3332:    Input Parameter:
3333: .  snes - the SNES context

3335:    Output Parameter:
3336: .  reason - negative value indicates diverged, positive value converged, see SNESConvergedReason or the
3337:             manual pages for the individual convergence tests for complete lists

3339:    Level: intermediate

3341:    Notes: Can only be called after the call the SNESSolve() is complete.

3343: .keywords: SNES, nonlinear, set, convergence, test

3345: .seealso: SNESSetConvergenceTest(), SNESConvergedReason
3346: @*/
3347: PetscErrorCode  SNESGetConvergedReason(SNES snes,SNESConvergedReason *reason)
3348: {
3352:   *reason = snes->reason;
3353:   return(0);
3354: }

3358: /*@
3359:    SNESSetConvergenceHistory - Sets the array used to hold the convergence history.

3361:    Logically Collective on SNES

3363:    Input Parameters:
3364: +  snes - iterative context obtained from SNESCreate()
3365: .  a   - array to hold history, this array will contain the function norms computed at each step
3366: .  its - integer array holds the number of linear iterations for each solve.
3367: .  na  - size of a and its
3368: -  reset - PETSC_TRUE indicates each new nonlinear solve resets the history counter to zero,
3369:            else it continues storing new values for new nonlinear solves after the old ones

3371:    Notes:
3372:    If 'a' and 'its' are NULL then space is allocated for the history. If 'na' PETSC_DECIDE or PETSC_DEFAULT then a
3373:    default array of length 10000 is allocated.

3375:    This routine is useful, e.g., when running a code for purposes
3376:    of accurate performance monitoring, when no I/O should be done
3377:    during the section of code that is being timed.

3379:    Level: intermediate

3381: .keywords: SNES, set, convergence, history

3383: .seealso: SNESGetConvergenceHistory()

3385: @*/
3386: PetscErrorCode  SNESSetConvergenceHistory(SNES snes,PetscReal a[],PetscInt its[],PetscInt na,PetscBool reset)
3387: {

3394:   if (!a) {
3395:     if (na == PETSC_DECIDE || na == PETSC_DEFAULT) na = 1000;
3396:     PetscMalloc(na*sizeof(PetscReal),&a);
3397:     PetscMalloc(na*sizeof(PetscInt),&its);

3399:     snes->conv_malloc = PETSC_TRUE;
3400:   }
3401:   snes->conv_hist       = a;
3402:   snes->conv_hist_its   = its;
3403:   snes->conv_hist_max   = na;
3404:   snes->conv_hist_len   = 0;
3405:   snes->conv_hist_reset = reset;
3406:   return(0);
3407: }

3409: #if defined(PETSC_HAVE_MATLAB_ENGINE)
3410: #include <engine.h>   /* MATLAB include file */
3411: #include <mex.h>      /* MATLAB include file */

3415: PETSC_EXTERN mxArray *SNESGetConvergenceHistoryMatlab(SNES snes)
3416: {
3417:   mxArray   *mat;
3418:   PetscInt  i;
3419:   PetscReal *ar;

3422:   mat = mxCreateDoubleMatrix(snes->conv_hist_len,1,mxREAL);
3423:   ar  = (PetscReal*) mxGetData(mat);
3424:   for (i=0; i<snes->conv_hist_len; i++) ar[i] = snes->conv_hist[i];
3425:   PetscFunctionReturn(mat);
3426: }
3427: #endif

3431: /*@C
3432:    SNESGetConvergenceHistory - Gets the array used to hold the convergence history.

3434:    Not Collective

3436:    Input Parameter:
3437: .  snes - iterative context obtained from SNESCreate()

3439:    Output Parameters:
3440: .  a   - array to hold history
3441: .  its - integer array holds the number of linear iterations (or
3442:          negative if not converged) for each solve.
3443: -  na  - size of a and its

3445:    Notes:
3446:     The calling sequence for this routine in Fortran is
3447: $   call SNESGetConvergenceHistory(SNES snes, integer na, integer ierr)

3449:    This routine is useful, e.g., when running a code for purposes
3450:    of accurate performance monitoring, when no I/O should be done
3451:    during the section of code that is being timed.

3453:    Level: intermediate

3455: .keywords: SNES, get, convergence, history

3457: .seealso: SNESSetConvergencHistory()

3459: @*/
3460: PetscErrorCode  SNESGetConvergenceHistory(SNES snes,PetscReal *a[],PetscInt *its[],PetscInt *na)
3461: {
3464:   if (a)   *a   = snes->conv_hist;
3465:   if (its) *its = snes->conv_hist_its;
3466:   if (na)  *na  = snes->conv_hist_len;
3467:   return(0);
3468: }

3472: /*@C
3473:   SNESSetUpdate - Sets the general-purpose update function called
3474:   at the beginning of every iteration of the nonlinear solve. Specifically
3475:   it is called just before the Jacobian is "evaluated".

3477:   Logically Collective on SNES

3479:   Input Parameters:
3480: . snes - The nonlinear solver context
3481: . func - The function

3483:   Calling sequence of func:
3484: . func (SNES snes, PetscInt step);

3486: . step - The current step of the iteration

3488:   Level: advanced

3490:   Note: This is NOT what one uses to update the ghost points before a function evaluation, that should be done at the beginning of your FormFunction()
3491:         This is not used by most users.

3493: .keywords: SNES, update

3495: .seealso SNESSetJacobian(), SNESSolve()
3496: @*/
3497: PetscErrorCode  SNESSetUpdate(SNES snes, PetscErrorCode (*func)(SNES, PetscInt))
3498: {
3501:   snes->ops->update = func;
3502:   return(0);
3503: }

3507: /*
3508:    SNESScaleStep_Private - Scales a step so that its length is less than the
3509:    positive parameter delta.

3511:     Input Parameters:
3512: +   snes - the SNES context
3513: .   y - approximate solution of linear system
3514: .   fnorm - 2-norm of current function
3515: -   delta - trust region size

3517:     Output Parameters:
3518: +   gpnorm - predicted function norm at the new point, assuming local
3519:     linearization.  The value is zero if the step lies within the trust
3520:     region, and exceeds zero otherwise.
3521: -   ynorm - 2-norm of the step

3523:     Note:
3524:     For non-trust region methods such as SNESNEWTONLS, the parameter delta
3525:     is set to be the maximum allowable step size.

3527: .keywords: SNES, nonlinear, scale, step
3528: */
3529: PetscErrorCode SNESScaleStep_Private(SNES snes,Vec y,PetscReal *fnorm,PetscReal *delta,PetscReal *gpnorm,PetscReal *ynorm)
3530: {
3531:   PetscReal      nrm;
3532:   PetscScalar    cnorm;


3540:   VecNorm(y,NORM_2,&nrm);
3541:   if (nrm > *delta) {
3542:     nrm     = *delta/nrm;
3543:     *gpnorm = (1.0 - nrm)*(*fnorm);
3544:     cnorm   = nrm;
3545:     VecScale(y,cnorm);
3546:     *ynorm  = *delta;
3547:   } else {
3548:     *gpnorm = 0.0;
3549:     *ynorm  = nrm;
3550:   }
3551:   return(0);
3552: }

3556: /*@C
3557:    SNESSolve - Solves a nonlinear system F(x) = b.
3558:    Call SNESSolve() after calling SNESCreate() and optional routines of the form SNESSetXXX().

3560:    Collective on SNES

3562:    Input Parameters:
3563: +  snes - the SNES context
3564: .  b - the constant part of the equation F(x) = b, or NULL to use zero.
3565: -  x - the solution vector.

3567:    Notes:
3568:    The user should initialize the vector,x, with the initial guess
3569:    for the nonlinear solve prior to calling SNESSolve.  In particular,
3570:    to employ an initial guess of zero, the user should explicitly set
3571:    this vector to zero by calling VecSet().

3573:    Level: beginner

3575: .keywords: SNES, nonlinear, solve

3577: .seealso: SNESCreate(), SNESDestroy(), SNESSetFunction(), SNESSetJacobian(), SNESSetGridSequence(), SNESGetSolution()
3578: @*/
3579: PetscErrorCode  SNESSolve(SNES snes,Vec b,Vec x)
3580: {
3581:   PetscErrorCode    ierr;
3582:   PetscBool         flg;
3583:   PetscViewer       viewer;
3584:   PetscInt          grid;
3585:   Vec               xcreated = NULL;
3586:   DM                dm;
3587:   PetscViewerFormat format;


3596:   if (!x) {
3597:     SNESGetDM(snes,&dm);
3598:     DMCreateGlobalVector(dm,&xcreated);
3599:     x    = xcreated;
3600:   }

3602:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"-snes_view_pre",&viewer,&format,&flg);
3603:   if (flg && !PetscPreLoadingOn) {
3604:     PetscViewerPushFormat(viewer,format);
3605:     SNESView(snes,viewer);
3606:     PetscViewerPopFormat(viewer);
3607:     PetscViewerDestroy(&viewer);
3608:   }

3610:   for (grid=0; grid<snes->gridsequence; grid++) {PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));}
3611:   for (grid=0; grid<snes->gridsequence+1; grid++) {

3613:     /* set solution vector */
3614:     if (!grid) {PetscObjectReference((PetscObject)x);}
3615:     VecDestroy(&snes->vec_sol);
3616:     snes->vec_sol = x;
3617:     SNESGetDM(snes,&dm);

3619:     /* set affine vector if provided */
3620:     if (b) { PetscObjectReference((PetscObject)b); }
3621:     VecDestroy(&snes->vec_rhs);
3622:     snes->vec_rhs = b;

3624:     SNESSetUp(snes);

3626:     if (!grid) {
3627:       if (snes->ops->computeinitialguess) {
3628:         (*snes->ops->computeinitialguess)(snes,snes->vec_sol,snes->initialguessP);
3629:       }
3630:     }

3632:     if (snes->conv_hist_reset) snes->conv_hist_len = 0;
3633:     snes->nfuncs = 0; snes->linear_its = 0; snes->numFailures = 0;

3635:     PetscLogEventBegin(SNES_Solve,snes,0,0,0);
3636:     (*snes->ops->solve)(snes);
3637:     PetscLogEventEnd(SNES_Solve,snes,0,0,0);
3638:     if (snes->domainerror) {
3639:       snes->reason      = SNES_DIVERGED_FUNCTION_DOMAIN;
3640:       snes->domainerror = PETSC_FALSE;
3641:     }
3642:     if (!snes->reason) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");

3644:     flg  = PETSC_FALSE;
3645:     PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_test_local_min",&flg,NULL);
3646:     if (flg && !PetscPreLoadingOn) { SNESTestLocalMin(snes); }
3647:     if (snes->printreason) {
3648:       PetscViewerASCIIAddTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),((PetscObject)snes)->tablevel);
3649:       if (snes->reason > 0) {
3650:         PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),"Nonlinear solve converged due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);
3651:       } else {
3652:         PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),"Nonlinear solve did not converge due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);
3653:       }
3654:       PetscViewerASCIISubtractTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),((PetscObject)snes)->tablevel);
3655:     }

3657:     if (snes->errorifnotconverged && snes->reason < 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged");
3658:     if (grid <  snes->gridsequence) {
3659:       DM  fine;
3660:       Vec xnew;
3661:       Mat interp;

3663:       DMRefine(snes->dm,PetscObjectComm((PetscObject)snes),&fine);
3664:       if (!fine) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"DMRefine() did not perform any refinement, cannot continue grid sequencing");
3665:       DMCreateInterpolation(snes->dm,fine,&interp,NULL);
3666:       DMCreateGlobalVector(fine,&xnew);
3667:       MatInterpolate(interp,x,xnew);
3668:       DMInterpolate(snes->dm,interp,fine);
3669:       MatDestroy(&interp);
3670:       x    = xnew;

3672:       SNESReset(snes);
3673:       SNESSetDM(snes,fine);
3674:       DMDestroy(&fine);
3675:       PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));
3676:     }
3677:   }
3678:   /* monitoring and viewing */
3679:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"-snes_view",&viewer,&format,&flg);
3680:   if (flg && !PetscPreLoadingOn) {
3681:     PetscViewerPushFormat(viewer,format);
3682:     SNESView(snes,viewer);
3683:     PetscViewerPopFormat(viewer);
3684:     PetscViewerDestroy(&viewer);
3685:   }
3686:   VecViewFromOptions(snes->vec_sol,((PetscObject)snes)->prefix,"-snes_view_solution");

3688:   VecDestroy(&xcreated);
3689:   PetscObjectAMSBlock((PetscObject)snes);
3690:   return(0);
3691: }

3693: /* --------- Internal routines for SNES Package --------- */

3697: /*@C
3698:    SNESSetType - Sets the method for the nonlinear solver.

3700:    Collective on SNES

3702:    Input Parameters:
3703: +  snes - the SNES context
3704: -  type - a known method

3706:    Options Database Key:
3707: .  -snes_type <type> - Sets the method; use -help for a list
3708:    of available methods (for instance, newtonls or newtontr)

3710:    Notes:
3711:    See "petsc/include/petscsnes.h" for available methods (for instance)
3712: +    SNESNEWTONLS - Newton's method with line search
3713:      (systems of nonlinear equations)
3714: .    SNESNEWTONTR - Newton's method with trust region
3715:      (systems of nonlinear equations)

3717:   Normally, it is best to use the SNESSetFromOptions() command and then
3718:   set the SNES solver type from the options database rather than by using
3719:   this routine.  Using the options database provides the user with
3720:   maximum flexibility in evaluating the many nonlinear solvers.
3721:   The SNESSetType() routine is provided for those situations where it
3722:   is necessary to set the nonlinear solver independently of the command
3723:   line or options database.  This might be the case, for example, when
3724:   the choice of solver changes during the execution of the program,
3725:   and the user's application is taking responsibility for choosing the
3726:   appropriate method.

3728:     Developer Notes: SNESRegister() adds a constructor for a new SNESType to SNESList, SNESSetType() locates
3729:     the constructor in that list and calls it to create the spexific object.

3731:   Level: intermediate

3733: .keywords: SNES, set, type

3735: .seealso: SNESType, SNESCreate(), SNESDestroy(), SNESGetType(), SNESSetFromOptions()

3737: @*/
3738: PetscErrorCode  SNESSetType(SNES snes,SNESType type)
3739: {
3740:   PetscErrorCode ierr,(*r)(SNES);
3741:   PetscBool      match;


3747:   PetscObjectTypeCompare((PetscObject)snes,type,&match);
3748:   if (match) return(0);

3750:    PetscFunctionListFind(SNESList,type,&r);
3751:   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested SNES type %s",type);
3752:   /* Destroy the previous private SNES context */
3753:   if (snes->ops->destroy) {
3754:     (*(snes)->ops->destroy)(snes);
3755:     snes->ops->destroy = NULL;
3756:   }
3757:   /* Reinitialize function pointers in SNESOps structure */
3758:   snes->ops->setup          = 0;
3759:   snes->ops->solve          = 0;
3760:   snes->ops->view           = 0;
3761:   snes->ops->setfromoptions = 0;
3762:   snes->ops->destroy        = 0;
3763:   /* Call the SNESCreate_XXX routine for this particular Nonlinear solver */
3764:   snes->setupcalled = PETSC_FALSE;

3766:   PetscObjectChangeTypeName((PetscObject)snes,type);
3767:   (*r)(snes);
3768:   return(0);
3769: }

3773: /*@C
3774:    SNESGetType - Gets the SNES method type and name (as a string).

3776:    Not Collective

3778:    Input Parameter:
3779: .  snes - nonlinear solver context

3781:    Output Parameter:
3782: .  type - SNES method (a character string)

3784:    Level: intermediate

3786: .keywords: SNES, nonlinear, get, type, name
3787: @*/
3788: PetscErrorCode  SNESGetType(SNES snes,SNESType *type)
3789: {
3793:   *type = ((PetscObject)snes)->type_name;
3794:   return(0);
3795: }

3799: /*@
3800:    SNESGetSolution - Returns the vector where the approximate solution is
3801:    stored. This is the fine grid solution when using SNESSetGridSequence().

3803:    Not Collective, but Vec is parallel if SNES is parallel

3805:    Input Parameter:
3806: .  snes - the SNES context

3808:    Output Parameter:
3809: .  x - the solution

3811:    Level: intermediate

3813: .keywords: SNES, nonlinear, get, solution

3815: .seealso:  SNESGetSolutionUpdate(), SNESGetFunction()
3816: @*/
3817: PetscErrorCode  SNESGetSolution(SNES snes,Vec *x)
3818: {
3822:   *x = snes->vec_sol;
3823:   return(0);
3824: }

3828: /*@
3829:    SNESGetSolutionUpdate - Returns the vector where the solution update is
3830:    stored.

3832:    Not Collective, but Vec is parallel if SNES is parallel

3834:    Input Parameter:
3835: .  snes - the SNES context

3837:    Output Parameter:
3838: .  x - the solution update

3840:    Level: advanced

3842: .keywords: SNES, nonlinear, get, solution, update

3844: .seealso: SNESGetSolution(), SNESGetFunction()
3845: @*/
3846: PetscErrorCode  SNESGetSolutionUpdate(SNES snes,Vec *x)
3847: {
3851:   *x = snes->vec_sol_update;
3852:   return(0);
3853: }

3857: /*@C
3858:    SNESGetFunction - Returns the vector where the function is stored.

3860:    Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet.

3862:    Input Parameter:
3863: .  snes - the SNES context

3865:    Output Parameter:
3866: +  r - the vector that is used to store residuals (or NULL if you don't want it)
3867: .  SNESFunction- the function (or NULL if you don't want it)
3868: -  ctx - the function context (or NULL if you don't want it)

3870:    Level: advanced

3872: .keywords: SNES, nonlinear, get, function

3874: .seealso: SNESSetFunction(), SNESGetSolution(), SNESFunction
3875: @*/
3876: PetscErrorCode  SNESGetFunction(SNES snes,Vec *r,PetscErrorCode (**SNESFunction)(SNES,Vec,Vec,void*),void **ctx)
3877: {
3879:   DM             dm;

3883:   if (r) {
3884:     if (!snes->vec_func) {
3885:       if (snes->vec_rhs) {
3886:         VecDuplicate(snes->vec_rhs,&snes->vec_func);
3887:       } else if (snes->vec_sol) {
3888:         VecDuplicate(snes->vec_sol,&snes->vec_func);
3889:       } else if (snes->dm) {
3890:         DMCreateGlobalVector(snes->dm,&snes->vec_func);
3891:       }
3892:     }
3893:     *r = snes->vec_func;
3894:   }
3895:   SNESGetDM(snes,&dm);
3896:   DMSNESGetFunction(dm,SNESFunction,ctx);
3897:   return(0);
3898: }

3900: /*@C
3901:    SNESGetGS - Returns the GS function and context.

3903:    Input Parameter:
3904: .  snes - the SNES context

3906:    Output Parameter:
3907: +  SNESGSFunction - the function (or NULL)
3908: -  ctx    - the function context (or NULL)

3910:    Level: advanced

3912: .keywords: SNES, nonlinear, get, function

3914: .seealso: SNESSetGS(), SNESGetFunction()
3915: @*/

3919: PetscErrorCode SNESGetGS (SNES snes, PetscErrorCode (**SNESGSFunction)(SNES, Vec, Vec, void*), void ** ctx)
3920: {
3922:   DM             dm;

3926:   SNESGetDM(snes,&dm);
3927:   DMSNESGetGS(dm,SNESGSFunction,ctx);
3928:   return(0);
3929: }

3933: /*@C
3934:    SNESSetOptionsPrefix - Sets the prefix used for searching for all
3935:    SNES options in the database.

3937:    Logically Collective on SNES

3939:    Input Parameter:
3940: +  snes - the SNES context
3941: -  prefix - the prefix to prepend to all option names

3943:    Notes:
3944:    A hyphen (-) must NOT be given at the beginning of the prefix name.
3945:    The first character of all runtime options is AUTOMATICALLY the hyphen.

3947:    Level: advanced

3949: .keywords: SNES, set, options, prefix, database

3951: .seealso: SNESSetFromOptions()
3952: @*/
3953: PetscErrorCode  SNESSetOptionsPrefix(SNES snes,const char prefix[])
3954: {

3959:   PetscObjectSetOptionsPrefix((PetscObject)snes,prefix);
3960:   if (!snes->ksp) {SNESGetKSP(snes,&snes->ksp);}
3961:   if (snes->linesearch) {
3962:     SNESGetLineSearch(snes,&snes->linesearch);
3963:     PetscObjectSetOptionsPrefix((PetscObject)snes->linesearch,prefix);
3964:   }
3965:   KSPSetOptionsPrefix(snes->ksp,prefix);
3966:   return(0);
3967: }

3971: /*@C
3972:    SNESAppendOptionsPrefix - Appends to the prefix used for searching for all
3973:    SNES options in the database.

3975:    Logically Collective on SNES

3977:    Input Parameters:
3978: +  snes - the SNES context
3979: -  prefix - the prefix to prepend to all option names

3981:    Notes:
3982:    A hyphen (-) must NOT be given at the beginning of the prefix name.
3983:    The first character of all runtime options is AUTOMATICALLY the hyphen.

3985:    Level: advanced

3987: .keywords: SNES, append, options, prefix, database

3989: .seealso: SNESGetOptionsPrefix()
3990: @*/
3991: PetscErrorCode  SNESAppendOptionsPrefix(SNES snes,const char prefix[])
3992: {

3997:   PetscObjectAppendOptionsPrefix((PetscObject)snes,prefix);
3998:   if (!snes->ksp) {SNESGetKSP(snes,&snes->ksp);}
3999:   if (snes->linesearch) {
4000:     SNESGetLineSearch(snes,&snes->linesearch);
4001:     PetscObjectAppendOptionsPrefix((PetscObject)snes->linesearch,prefix);
4002:   }
4003:   KSPAppendOptionsPrefix(snes->ksp,prefix);
4004:   return(0);
4005: }

4009: /*@C
4010:    SNESGetOptionsPrefix - Sets the prefix used for searching for all
4011:    SNES options in the database.

4013:    Not Collective

4015:    Input Parameter:
4016: .  snes - the SNES context

4018:    Output Parameter:
4019: .  prefix - pointer to the prefix string used

4021:    Notes: On the fortran side, the user should pass in a string 'prefix' of
4022:    sufficient length to hold the prefix.

4024:    Level: advanced

4026: .keywords: SNES, get, options, prefix, database

4028: .seealso: SNESAppendOptionsPrefix()
4029: @*/
4030: PetscErrorCode  SNESGetOptionsPrefix(SNES snes,const char *prefix[])
4031: {

4036:   PetscObjectGetOptionsPrefix((PetscObject)snes,prefix);
4037:   return(0);
4038: }


4043: /*@C
4044:   SNESRegister - Adds a method to the nonlinear solver package.

4046:    Not collective

4048:    Input Parameters:
4049: +  name_solver - name of a new user-defined solver
4050: -  routine_create - routine to create method context

4052:    Notes:
4053:    SNESRegister() may be called multiple times to add several user-defined solvers.

4055:    Sample usage:
4056: .vb
4057:    SNESRegister("my_solver",MySolverCreate);
4058: .ve

4060:    Then, your solver can be chosen with the procedural interface via
4061: $     SNESSetType(snes,"my_solver")
4062:    or at runtime via the option
4063: $     -snes_type my_solver

4065:    Level: advanced

4067:     Note: If your function is not being put into a shared library then use SNESRegister() instead

4069: .keywords: SNES, nonlinear, register

4071: .seealso: SNESRegisterAll(), SNESRegisterDestroy()

4073:   Level: advanced
4074: @*/
4075: PetscErrorCode  SNESRegister(const char sname[],PetscErrorCode (*function)(SNES))
4076: {

4080:   PetscFunctionListAdd(&SNESList,sname,function);
4081:   return(0);
4082: }

4086: PetscErrorCode  SNESTestLocalMin(SNES snes)
4087: {
4089:   PetscInt       N,i,j;
4090:   Vec            u,uh,fh;
4091:   PetscScalar    value;
4092:   PetscReal      norm;

4095:   SNESGetSolution(snes,&u);
4096:   VecDuplicate(u,&uh);
4097:   VecDuplicate(u,&fh);

4099:   /* currently only works for sequential */
4100:   PetscPrintf(PETSC_COMM_WORLD,"Testing FormFunction() for local min\n");
4101:   VecGetSize(u,&N);
4102:   for (i=0; i<N; i++) {
4103:     VecCopy(u,uh);
4104:     PetscPrintf(PETSC_COMM_WORLD,"i = %D\n",i);
4105:     for (j=-10; j<11; j++) {
4106:       value = PetscSign(j)*exp(PetscAbs(j)-10.0);
4107:       VecSetValue(uh,i,value,ADD_VALUES);
4108:       SNESComputeFunction(snes,uh,fh);
4109:       VecNorm(fh,NORM_2,&norm);
4110:       PetscPrintf(PETSC_COMM_WORLD,"       j norm %D %18.16e\n",j,norm);
4111:       value = -value;
4112:       VecSetValue(uh,i,value,ADD_VALUES);
4113:     }
4114:   }
4115:   VecDestroy(&uh);
4116:   VecDestroy(&fh);
4117:   return(0);
4118: }

4122: /*@
4123:    SNESKSPSetUseEW - Sets SNES use Eisenstat-Walker method for
4124:    computing relative tolerance for linear solvers within an inexact
4125:    Newton method.

4127:    Logically Collective on SNES

4129:    Input Parameters:
4130: +  snes - SNES context
4131: -  flag - PETSC_TRUE or PETSC_FALSE

4133:     Options Database:
4134: +  -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence
4135: .  -snes_ksp_ew_version ver - version of  Eisenstat-Walker method
4136: .  -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0
4137: .  -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax
4138: .  -snes_ksp_ew_gamma <gamma> - Sets gamma
4139: .  -snes_ksp_ew_alpha <alpha> - Sets alpha
4140: .  -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2
4141: -  -snes_ksp_ew_threshold <threshold> - Sets threshold

4143:    Notes:
4144:    Currently, the default is to use a constant relative tolerance for
4145:    the inner linear solvers.  Alternatively, one can use the
4146:    Eisenstat-Walker method, where the relative convergence tolerance
4147:    is reset at each Newton iteration according progress of the nonlinear
4148:    solver.

4150:    Level: advanced

4152:    Reference:
4153:    S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an
4154:    inexact Newton method", SISC 17 (1), pp.16-32, 1996.

4156: .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton

4158: .seealso: SNESKSPGetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW()
4159: @*/
4160: PetscErrorCode  SNESKSPSetUseEW(SNES snes,PetscBool flag)
4161: {
4165:   snes->ksp_ewconv = flag;
4166:   return(0);
4167: }

4171: /*@
4172:    SNESKSPGetUseEW - Gets if SNES is using Eisenstat-Walker method
4173:    for computing relative tolerance for linear solvers within an
4174:    inexact Newton method.

4176:    Not Collective

4178:    Input Parameter:
4179: .  snes - SNES context

4181:    Output Parameter:
4182: .  flag - PETSC_TRUE or PETSC_FALSE

4184:    Notes:
4185:    Currently, the default is to use a constant relative tolerance for
4186:    the inner linear solvers.  Alternatively, one can use the
4187:    Eisenstat-Walker method, where the relative convergence tolerance
4188:    is reset at each Newton iteration according progress of the nonlinear
4189:    solver.

4191:    Level: advanced

4193:    Reference:
4194:    S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an
4195:    inexact Newton method", SISC 17 (1), pp.16-32, 1996.

4197: .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton

4199: .seealso: SNESKSPSetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW()
4200: @*/
4201: PetscErrorCode  SNESKSPGetUseEW(SNES snes, PetscBool  *flag)
4202: {
4206:   *flag = snes->ksp_ewconv;
4207:   return(0);
4208: }

4212: /*@
4213:    SNESKSPSetParametersEW - Sets parameters for Eisenstat-Walker
4214:    convergence criteria for the linear solvers within an inexact
4215:    Newton method.

4217:    Logically Collective on SNES

4219:    Input Parameters:
4220: +    snes - SNES context
4221: .    version - version 1, 2 (default is 2) or 3
4222: .    rtol_0 - initial relative tolerance (0 <= rtol_0 < 1)
4223: .    rtol_max - maximum relative tolerance (0 <= rtol_max < 1)
4224: .    gamma - multiplicative factor for version 2 rtol computation
4225:              (0 <= gamma2 <= 1)
4226: .    alpha - power for version 2 rtol computation (1 < alpha <= 2)
4227: .    alpha2 - power for safeguard
4228: -    threshold - threshold for imposing safeguard (0 < threshold < 1)

4230:    Note:
4231:    Version 3 was contributed by Luis Chacon, June 2006.

4233:    Use PETSC_DEFAULT to retain the default for any of the parameters.

4235:    Level: advanced

4237:    Reference:
4238:    S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an
4239:    inexact Newton method", Utah State University Math. Stat. Dept. Res.
4240:    Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput.

4242: .keywords: SNES, KSP, Eisenstat, Walker, set, parameters

4244: .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPGetParametersEW()
4245: @*/
4246: PetscErrorCode  SNESKSPSetParametersEW(SNES snes,PetscInt version,PetscReal rtol_0,PetscReal rtol_max,PetscReal gamma,PetscReal alpha,PetscReal alpha2,PetscReal threshold)
4247: {
4248:   SNESKSPEW *kctx;

4252:   kctx = (SNESKSPEW*)snes->kspconvctx;
4253:   if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing");

4262:   if (version != PETSC_DEFAULT)   kctx->version   = version;
4263:   if (rtol_0 != PETSC_DEFAULT)    kctx->rtol_0    = rtol_0;
4264:   if (rtol_max != PETSC_DEFAULT)  kctx->rtol_max  = rtol_max;
4265:   if (gamma != PETSC_DEFAULT)     kctx->gamma     = gamma;
4266:   if (alpha != PETSC_DEFAULT)     kctx->alpha     = alpha;
4267:   if (alpha2 != PETSC_DEFAULT)    kctx->alpha2    = alpha2;
4268:   if (threshold != PETSC_DEFAULT) kctx->threshold = threshold;

4270:   if (kctx->version < 1 || kctx->version > 3) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 and 3 are supported: %D",kctx->version);
4271:   if (kctx->rtol_0 < 0.0 || kctx->rtol_0 >= 1.0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= rtol_0 < 1.0: %G",kctx->rtol_0);
4272:   if (kctx->rtol_max < 0.0 || kctx->rtol_max >= 1.0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= rtol_max (%G) < 1.0\n",kctx->rtol_max);
4273:   if (kctx->gamma < 0.0 || kctx->gamma > 1.0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= gamma (%G) <= 1.0\n",kctx->gamma);
4274:   if (kctx->alpha <= 1.0 || kctx->alpha > 2.0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"1.0 < alpha (%G) <= 2.0\n",kctx->alpha);
4275:   if (kctx->threshold <= 0.0 || kctx->threshold >= 1.0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 < threshold (%G) < 1.0\n",kctx->threshold);
4276:   return(0);
4277: }

4281: /*@
4282:    SNESKSPGetParametersEW - Gets parameters for Eisenstat-Walker
4283:    convergence criteria for the linear solvers within an inexact
4284:    Newton method.

4286:    Not Collective

4288:    Input Parameters:
4289:      snes - SNES context

4291:    Output Parameters:
4292: +    version - version 1, 2 (default is 2) or 3
4293: .    rtol_0 - initial relative tolerance (0 <= rtol_0 < 1)
4294: .    rtol_max - maximum relative tolerance (0 <= rtol_max < 1)
4295: .    gamma - multiplicative factor for version 2 rtol computation (0 <= gamma2 <= 1)
4296: .    alpha - power for version 2 rtol computation (1 < alpha <= 2)
4297: .    alpha2 - power for safeguard
4298: -    threshold - threshold for imposing safeguard (0 < threshold < 1)

4300:    Level: advanced

4302: .keywords: SNES, KSP, Eisenstat, Walker, get, parameters

4304: .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPSetParametersEW()
4305: @*/
4306: PetscErrorCode  SNESKSPGetParametersEW(SNES snes,PetscInt *version,PetscReal *rtol_0,PetscReal *rtol_max,PetscReal *gamma,PetscReal *alpha,PetscReal *alpha2,PetscReal *threshold)
4307: {
4308:   SNESKSPEW *kctx;

4312:   kctx = (SNESKSPEW*)snes->kspconvctx;
4313:   if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing");
4314:   if (version)   *version   = kctx->version;
4315:   if (rtol_0)    *rtol_0    = kctx->rtol_0;
4316:   if (rtol_max)  *rtol_max  = kctx->rtol_max;
4317:   if (gamma)     *gamma     = kctx->gamma;
4318:   if (alpha)     *alpha     = kctx->alpha;
4319:   if (alpha2)    *alpha2    = kctx->alpha2;
4320:   if (threshold) *threshold = kctx->threshold;
4321:   return(0);
4322: }

4326:  PetscErrorCode SNESKSPEW_PreSolve(KSP ksp, Vec b, Vec x, SNES snes)
4327: {
4329:   SNESKSPEW      *kctx = (SNESKSPEW*)snes->kspconvctx;
4330:   PetscReal      rtol  = PETSC_DEFAULT,stol;

4333:   if (!snes->ksp_ewconv) return(0);
4334:   if (!snes->iter) rtol = kctx->rtol_0; /* first time in, so use the original user rtol */
4335:   else {
4336:     if (kctx->version == 1) {
4337:       rtol = (snes->norm - kctx->lresid_last)/kctx->norm_last;
4338:       if (rtol < 0.0) rtol = -rtol;
4339:       stol = PetscPowReal(kctx->rtol_last,kctx->alpha2);
4340:       if (stol > kctx->threshold) rtol = PetscMax(rtol,stol);
4341:     } else if (kctx->version == 2) {
4342:       rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha);
4343:       stol = kctx->gamma * PetscPowReal(kctx->rtol_last,kctx->alpha);
4344:       if (stol > kctx->threshold) rtol = PetscMax(rtol,stol);
4345:     } else if (kctx->version == 3) { /* contributed by Luis Chacon, June 2006. */
4346:       rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha);
4347:       /* safeguard: avoid sharp decrease of rtol */
4348:       stol = kctx->gamma*PetscPowReal(kctx->rtol_last,kctx->alpha);
4349:       stol = PetscMax(rtol,stol);
4350:       rtol = PetscMin(kctx->rtol_0,stol);
4351:       /* safeguard: avoid oversolving */
4352:       stol = kctx->gamma*(snes->ttol)/snes->norm;
4353:       stol = PetscMax(rtol,stol);
4354:       rtol = PetscMin(kctx->rtol_0,stol);
4355:     } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 or 3 are supported: %D",kctx->version);
4356:   }
4357:   /* safeguard: avoid rtol greater than one */
4358:   rtol = PetscMin(rtol,kctx->rtol_max);
4359:   KSPSetTolerances(ksp,rtol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
4360:   PetscInfo3(snes,"iter %D, Eisenstat-Walker (version %D) KSP rtol=%G\n",snes->iter,kctx->version,rtol);
4361:   return(0);
4362: }

4366: PetscErrorCode SNESKSPEW_PostSolve(KSP ksp, Vec b, Vec x, SNES snes)
4367: {
4369:   SNESKSPEW      *kctx = (SNESKSPEW*)snes->kspconvctx;
4370:   PCSide         pcside;
4371:   Vec            lres;

4374:   if (!snes->ksp_ewconv) return(0);
4375:   KSPGetTolerances(ksp,&kctx->rtol_last,0,0,0);
4376:   SNESGetFunctionNorm(snes,&kctx->norm_last);
4377:   if (kctx->version == 1) {
4378:     KSPGetPCSide(ksp,&pcside);
4379:     if (pcside == PC_RIGHT) { /* XXX Should we also test KSP_UNPRECONDITIONED_NORM ? */
4380:       /* KSP residual is true linear residual */
4381:       KSPGetResidualNorm(ksp,&kctx->lresid_last);
4382:     } else {
4383:       /* KSP residual is preconditioned residual */
4384:       /* compute true linear residual norm */
4385:       VecDuplicate(b,&lres);
4386:       MatMult(snes->jacobian,x,lres);
4387:       VecAYPX(lres,-1.0,b);
4388:       VecNorm(lres,NORM_2,&kctx->lresid_last);
4389:       VecDestroy(&lres);
4390:     }
4391:   }
4392:   return(0);
4393: }

4397: /*@
4398:    SNESGetKSP - Returns the KSP context for a SNES solver.

4400:    Not Collective, but if SNES object is parallel, then KSP object is parallel

4402:    Input Parameter:
4403: .  snes - the SNES context

4405:    Output Parameter:
4406: .  ksp - the KSP context

4408:    Notes:
4409:    The user can then directly manipulate the KSP context to set various
4410:    options, etc.  Likewise, the user can then extract and manipulate the
4411:    PC contexts as well.

4413:    Level: beginner

4415: .keywords: SNES, nonlinear, get, KSP, context

4417: .seealso: KSPGetPC(), SNESCreate(), KSPCreate(), SNESSetKSP()
4418: @*/
4419: PetscErrorCode  SNESGetKSP(SNES snes,KSP *ksp)
4420: {


4427:   if (!snes->ksp) {
4428:     KSPCreate(PetscObjectComm((PetscObject)snes),&snes->ksp);
4429:     PetscObjectIncrementTabLevel((PetscObject)snes->ksp,(PetscObject)snes,1);
4430:     PetscLogObjectParent(snes,snes->ksp);

4432:     KSPSetPreSolve(snes->ksp,(PetscErrorCode (*)(KSP,Vec,Vec,void*))SNESKSPEW_PreSolve,snes);
4433:     KSPSetPostSolve(snes->ksp,(PetscErrorCode (*)(KSP,Vec,Vec,void*))SNESKSPEW_PostSolve,snes);
4434:   }
4435:   *ksp = snes->ksp;
4436:   return(0);
4437: }


4440: #include <petsc-private/dmimpl.h>
4443: /*@
4444:    SNESSetDM - Sets the DM that may be used by some preconditioners

4446:    Logically Collective on SNES

4448:    Input Parameters:
4449: +  snes - the preconditioner context
4450: -  dm - the dm

4452:    Level: intermediate

4454: .seealso: SNESGetDM(), KSPSetDM(), KSPGetDM()
4455: @*/
4456: PetscErrorCode  SNESSetDM(SNES snes,DM dm)
4457: {
4459:   KSP            ksp;
4460:   DMSNES         sdm;

4464:   if (dm) {PetscObjectReference((PetscObject)dm);}
4465:   if (snes->dm) {               /* Move the DMSNES context over to the new DM unless the new DM already has one */
4466:     if (snes->dm->dmsnes && snes->dmAuto && !dm->dmsnes) {
4467:       DMCopyDMSNES(snes->dm,dm);
4468:       DMGetDMSNES(snes->dm,&sdm);
4469:       if (sdm->originaldm == snes->dm) sdm->originaldm = dm; /* Grant write privileges to the replacement DM */
4470:     }
4471:     DMDestroy(&snes->dm);
4472:   }
4473:   snes->dm     = dm;
4474:   snes->dmAuto = PETSC_FALSE;

4476:   SNESGetKSP(snes,&ksp);
4477:   KSPSetDM(ksp,dm);
4478:   KSPSetDMActive(ksp,PETSC_FALSE);
4479:   if (snes->pc) {
4480:     SNESSetDM(snes->pc, snes->dm);
4481:     SNESSetPCSide(snes,snes->pcside);
4482:   }
4483:   return(0);
4484: }

4488: /*@
4489:    SNESGetDM - Gets the DM that may be used by some preconditioners

4491:    Not Collective but DM obtained is parallel on SNES

4493:    Input Parameter:
4494: . snes - the preconditioner context

4496:    Output Parameter:
4497: .  dm - the dm

4499:    Level: intermediate

4501: .seealso: SNESSetDM(), KSPSetDM(), KSPGetDM()
4502: @*/
4503: PetscErrorCode  SNESGetDM(SNES snes,DM *dm)
4504: {

4509:   if (!snes->dm) {
4510:     DMShellCreate(PetscObjectComm((PetscObject)snes),&snes->dm);
4511:     snes->dmAuto = PETSC_TRUE;
4512:   }
4513:   *dm = snes->dm;
4514:   return(0);
4515: }

4519: /*@
4520:   SNESSetPC - Sets the nonlinear preconditioner to be used.

4522:   Collective on SNES

4524:   Input Parameters:
4525: + snes - iterative context obtained from SNESCreate()
4526: - pc   - the preconditioner object

4528:   Notes:
4529:   Use SNESGetPC() to retrieve the preconditioner context (for example,
4530:   to configure it using the API).

4532:   Level: developer

4534: .keywords: SNES, set, precondition
4535: .seealso: SNESGetPC()
4536: @*/
4537: PetscErrorCode SNESSetPC(SNES snes, SNES pc)
4538: {

4545:   PetscObjectReference((PetscObject) pc);
4546:   SNESDestroy(&snes->pc);
4547:   snes->pc = pc;
4548:   PetscLogObjectParent(snes, snes->pc);
4549:   return(0);
4550: }

4554: /*@
4555:   SNESGetPC - Returns a pointer to the nonlinear preconditioning context set with SNESSetPC().

4557:   Not Collective

4559:   Input Parameter:
4560: . snes - iterative context obtained from SNESCreate()

4562:   Output Parameter:
4563: . pc - preconditioner context

4565:   Level: developer

4567: .keywords: SNES, get, preconditioner
4568: .seealso: SNESSetPC()
4569: @*/
4570: PetscErrorCode SNESGetPC(SNES snes, SNES *pc)
4571: {
4573:   const char     *optionsprefix;

4578:   if (!snes->pc) {
4579:     SNESCreate(PetscObjectComm((PetscObject)snes),&snes->pc);
4580:     PetscObjectIncrementTabLevel((PetscObject)snes->pc,(PetscObject)snes,1);
4581:     PetscLogObjectParent(snes,snes->pc);
4582:     SNESGetOptionsPrefix(snes,&optionsprefix);
4583:     SNESSetOptionsPrefix(snes->pc,optionsprefix);
4584:     SNESAppendOptionsPrefix(snes->pc,"npc_");
4585:   }
4586:   *pc = snes->pc;
4587:   return(0);
4588: }

4592: /*@
4593:     SNESSetPCSide - Sets the preconditioning side.

4595:     Logically Collective on SNES

4597:     Input Parameter:
4598: .   snes - iterative context obtained from SNESCreate()

4600:     Output Parameter:
4601: .   side - the preconditioning side, where side is one of
4602: .vb
4603:       PC_LEFT - left preconditioning (default)
4604:       PC_RIGHT - right preconditioning
4605: .ve

4607:     Options Database Keys:
4608: .   -snes_pc_side <right,left>

4610:     Level: intermediate

4612: .keywords: SNES, set, right, left, side, preconditioner, flag

4614: .seealso: SNESGetPCSide(), KSPSetPCSide()
4615: @*/
4616: PetscErrorCode  SNESSetPCSide(SNES snes,PCSide side)
4617: {
4621:   snes->pcside = side;
4622:   return(0);
4623: }

4627: /*@
4628:     SNESGetPCSide - Gets the preconditioning side.

4630:     Not Collective

4632:     Input Parameter:
4633: .   snes - iterative context obtained from SNESCreate()

4635:     Output Parameter:
4636: .   side - the preconditioning side, where side is one of
4637: .vb
4638:       PC_LEFT - left preconditioning (default)
4639:       PC_RIGHT - right preconditioning
4640: .ve

4642:     Level: intermediate

4644: .keywords: SNES, get, right, left, side, preconditioner, flag

4646: .seealso: SNESSetPCSide(), KSPGetPCSide()
4647: @*/
4648: PetscErrorCode  SNESGetPCSide(SNES snes,PCSide *side)
4649: {
4653:   *side = snes->pcside;
4654:   return(0);
4655: }

4659: /*@
4660:   SNESSetLineSearch - Sets the linesearch on the SNES instance.

4662:   Collective on SNES

4664:   Input Parameters:
4665: + snes - iterative context obtained from SNESCreate()
4666: - linesearch   - the linesearch object

4668:   Notes:
4669:   Use SNESGetLineSearch() to retrieve the preconditioner context (for example,
4670:   to configure it using the API).

4672:   Level: developer

4674: .keywords: SNES, set, linesearch
4675: .seealso: SNESGetLineSearch()
4676: @*/
4677: PetscErrorCode SNESSetLineSearch(SNES snes, SNESLineSearch linesearch)
4678: {

4685:   PetscObjectReference((PetscObject) linesearch);
4686:   SNESLineSearchDestroy(&snes->linesearch);

4688:   snes->linesearch = linesearch;

4690:   PetscLogObjectParent(snes, snes->linesearch);
4691:   return(0);
4692: }

4696: /*@
4697:   SNESGetLineSearch - Returns a pointer to the line search context set with SNESSetLineSearch()
4698:   or creates a default line search instance associated with the SNES and returns it.

4700:   Not Collective

4702:   Input Parameter:
4703: . snes - iterative context obtained from SNESCreate()

4705:   Output Parameter:
4706: . linesearch - linesearch context

4708:   Level: developer

4710: .keywords: SNES, get, linesearch
4711: .seealso: SNESSetLineSearch()
4712: @*/
4713: PetscErrorCode SNESGetLineSearch(SNES snes, SNESLineSearch *linesearch)
4714: {
4716:   const char     *optionsprefix;

4721:   if (!snes->linesearch) {
4722:     SNESGetOptionsPrefix(snes, &optionsprefix);
4723:     SNESLineSearchCreate(PetscObjectComm((PetscObject)snes), &snes->linesearch);
4724:     SNESLineSearchSetSNES(snes->linesearch, snes);
4725:     SNESLineSearchAppendOptionsPrefix(snes->linesearch, optionsprefix);
4726:     PetscObjectIncrementTabLevel((PetscObject) snes->linesearch, (PetscObject) snes, 1);
4727:     PetscLogObjectParent(snes, snes->linesearch);
4728:   }
4729:   *linesearch = snes->linesearch;
4730:   return(0);
4731: }

4733: #if defined(PETSC_HAVE_MATLAB_ENGINE)
4734: #include <mex.h>

4736: typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext;

4740: /*
4741:    SNESComputeFunction_Matlab - Calls the function that has been set with SNESSetFunctionMatlab().

4743:    Collective on SNES

4745:    Input Parameters:
4746: +  snes - the SNES context
4747: -  x - input vector

4749:    Output Parameter:
4750: .  y - function vector, as set by SNESSetFunction()

4752:    Notes:
4753:    SNESComputeFunction() is typically used within nonlinear solvers
4754:    implementations, so most users would not generally call this routine
4755:    themselves.

4757:    Level: developer

4759: .keywords: SNES, nonlinear, compute, function

4761: .seealso: SNESSetFunction(), SNESGetFunction()
4762: */
4763: PetscErrorCode  SNESComputeFunction_Matlab(SNES snes,Vec x,Vec y, void *ctx)
4764: {
4765:   PetscErrorCode    ierr;
4766:   SNESMatlabContext *sctx = (SNESMatlabContext*)ctx;
4767:   int               nlhs  = 1,nrhs = 5;
4768:   mxArray           *plhs[1],*prhs[5];
4769:   long long int     lx = 0,ly = 0,ls = 0;


4778:   /* call Matlab function in ctx with arguments x and y */

4780:   PetscMemcpy(&ls,&snes,sizeof(snes));
4781:   PetscMemcpy(&lx,&x,sizeof(x));
4782:   PetscMemcpy(&ly,&y,sizeof(x));
4783:   prhs[0] = mxCreateDoubleScalar((double)ls);
4784:   prhs[1] = mxCreateDoubleScalar((double)lx);
4785:   prhs[2] = mxCreateDoubleScalar((double)ly);
4786:   prhs[3] = mxCreateString(sctx->funcname);
4787:   prhs[4] = sctx->ctx;
4788:   mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeFunctionInternal");
4789:   mxGetScalar(plhs[0]);
4790:   mxDestroyArray(prhs[0]);
4791:   mxDestroyArray(prhs[1]);
4792:   mxDestroyArray(prhs[2]);
4793:   mxDestroyArray(prhs[3]);
4794:   mxDestroyArray(plhs[0]);
4795:   return(0);
4796: }

4800: /*
4801:    SNESSetFunctionMatlab - Sets the function evaluation routine and function
4802:    vector for use by the SNES routines in solving systems of nonlinear
4803:    equations from MATLAB. Here the function is a string containing the name of a MATLAB function

4805:    Logically Collective on SNES

4807:    Input Parameters:
4808: +  snes - the SNES context
4809: .  r - vector to store function value
4810: -  SNESFunction - function evaluation routine

4812:    Notes:
4813:    The Newton-like methods typically solve linear systems of the form
4814: $      f'(x) x = -f(x),
4815:    where f'(x) denotes the Jacobian matrix and f(x) is the function.

4817:    Level: beginner

4819:    Developer Note:  This bleeds the allocated memory SNESMatlabContext *sctx;

4821: .keywords: SNES, nonlinear, set, function

4823: .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction()
4824: */
4825: PetscErrorCode  SNESSetFunctionMatlab(SNES snes,Vec r,const char *SNESFunction,mxArray *ctx)
4826: {
4827:   PetscErrorCode    ierr;
4828:   SNESMatlabContext *sctx;

4831:   /* currently sctx is memory bleed */
4832:   PetscMalloc(sizeof(SNESMatlabContext),&sctx);
4833:   PetscStrallocpy(SNESFunction,&sctx->funcname);
4834:   /*
4835:      This should work, but it doesn't
4836:   sctx->ctx = ctx;
4837:   mexMakeArrayPersistent(sctx->ctx);
4838:   */
4839:   sctx->ctx = mxDuplicateArray(ctx);
4840:   SNESSetFunction(snes,r,SNESComputeFunction_Matlab,sctx);
4841:   return(0);
4842: }

4846: /*
4847:    SNESComputeJacobian_Matlab - Calls the function that has been set with SNESSetJacobianMatlab().

4849:    Collective on SNES

4851:    Input Parameters:
4852: +  snes - the SNES context
4853: .  x - input vector
4854: .  A, B - the matrices
4855: -  ctx - user context

4857:    Output Parameter:
4858: .  flag - structure of the matrix

4860:    Level: developer

4862: .keywords: SNES, nonlinear, compute, function

4864: .seealso: SNESSetFunction(), SNESGetFunction()
4865: @*/
4866: PetscErrorCode  SNESComputeJacobian_Matlab(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *flag, void *ctx)
4867: {
4868:   PetscErrorCode    ierr;
4869:   SNESMatlabContext *sctx = (SNESMatlabContext*)ctx;
4870:   int               nlhs  = 2,nrhs = 6;
4871:   mxArray           *plhs[2],*prhs[6];
4872:   long long int     lx = 0,lA = 0,ls = 0, lB = 0;


4878:   /* call Matlab function in ctx with arguments x and y */

4880:   PetscMemcpy(&ls,&snes,sizeof(snes));
4881:   PetscMemcpy(&lx,&x,sizeof(x));
4882:   PetscMemcpy(&lA,A,sizeof(x));
4883:   PetscMemcpy(&lB,B,sizeof(x));
4884:   prhs[0] = mxCreateDoubleScalar((double)ls);
4885:   prhs[1] = mxCreateDoubleScalar((double)lx);
4886:   prhs[2] = mxCreateDoubleScalar((double)lA);
4887:   prhs[3] = mxCreateDoubleScalar((double)lB);
4888:   prhs[4] = mxCreateString(sctx->funcname);
4889:   prhs[5] = sctx->ctx;
4890:   mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeJacobianInternal");
4891:   mxGetScalar(plhs[0]);
4892:   *flag   = (MatStructure) mxGetScalar(plhs[1]);
4893:   mxDestroyArray(prhs[0]);
4894:   mxDestroyArray(prhs[1]);
4895:   mxDestroyArray(prhs[2]);
4896:   mxDestroyArray(prhs[3]);
4897:   mxDestroyArray(prhs[4]);
4898:   mxDestroyArray(plhs[0]);
4899:   mxDestroyArray(plhs[1]);
4900:   return(0);
4901: }

4905: /*
4906:    SNESSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices
4907:    vector for use by the SNES routines in solving systems of nonlinear
4908:    equations from MATLAB. Here the function is a string containing the name of a MATLAB function

4910:    Logically Collective on SNES

4912:    Input Parameters:
4913: +  snes - the SNES context
4914: .  A,B - Jacobian matrices
4915: .  SNESJacobianFunction - function evaluation routine
4916: -  ctx - user context

4918:    Level: developer

4920:    Developer Note:  This bleeds the allocated memory SNESMatlabContext *sctx;

4922: .keywords: SNES, nonlinear, set, function

4924: .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction(), SNESJacobianFunction
4925: */
4926: PetscErrorCode  SNESSetJacobianMatlab(SNES snes,Mat A,Mat B,const char *SNESJacobianFunction,mxArray *ctx)
4927: {
4928:   PetscErrorCode    ierr;
4929:   SNESMatlabContext *sctx;

4932:   /* currently sctx is memory bleed */
4933:   PetscMalloc(sizeof(SNESMatlabContext),&sctx);
4934:   PetscStrallocpy(SNESJacobianFunction,&sctx->funcname);
4935:   /*
4936:      This should work, but it doesn't
4937:   sctx->ctx = ctx;
4938:   mexMakeArrayPersistent(sctx->ctx);
4939:   */
4940:   sctx->ctx = mxDuplicateArray(ctx);
4941:   SNESSetJacobian(snes,A,B,SNESComputeJacobian_Matlab,sctx);
4942:   return(0);
4943: }

4947: /*
4948:    SNESMonitor_Matlab - Calls the function that has been set with SNESMonitorSetMatlab().

4950:    Collective on SNES

4952: .seealso: SNESSetFunction(), SNESGetFunction()
4953: @*/
4954: PetscErrorCode  SNESMonitor_Matlab(SNES snes,PetscInt it, PetscReal fnorm, void *ctx)
4955: {
4956:   PetscErrorCode    ierr;
4957:   SNESMatlabContext *sctx = (SNESMatlabContext*)ctx;
4958:   int               nlhs  = 1,nrhs = 6;
4959:   mxArray           *plhs[1],*prhs[6];
4960:   long long int     lx = 0,ls = 0;
4961:   Vec               x  = snes->vec_sol;


4966:   PetscMemcpy(&ls,&snes,sizeof(snes));
4967:   PetscMemcpy(&lx,&x,sizeof(x));
4968:   prhs[0] = mxCreateDoubleScalar((double)ls);
4969:   prhs[1] = mxCreateDoubleScalar((double)it);
4970:   prhs[2] = mxCreateDoubleScalar((double)fnorm);
4971:   prhs[3] = mxCreateDoubleScalar((double)lx);
4972:   prhs[4] = mxCreateString(sctx->funcname);
4973:   prhs[5] = sctx->ctx;
4974:   mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESMonitorInternal");
4975:   mxGetScalar(plhs[0]);
4976:   mxDestroyArray(prhs[0]);
4977:   mxDestroyArray(prhs[1]);
4978:   mxDestroyArray(prhs[2]);
4979:   mxDestroyArray(prhs[3]);
4980:   mxDestroyArray(prhs[4]);
4981:   mxDestroyArray(plhs[0]);
4982:   return(0);
4983: }

4987: /*
4988:    SNESMonitorSetMatlab - Sets the monitor function from MATLAB

4990:    Level: developer

4992:    Developer Note:  This bleeds the allocated memory SNESMatlabContext *sctx;

4994: .keywords: SNES, nonlinear, set, function

4996: .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction()
4997: */
4998: PetscErrorCode  SNESMonitorSetMatlab(SNES snes,const char *SNESMonitorFunction,mxArray *ctx)
4999: {
5000:   PetscErrorCode    ierr;
5001:   SNESMatlabContext *sctx;

5004:   /* currently sctx is memory bleed */
5005:   PetscMalloc(sizeof(SNESMatlabContext),&sctx);
5006:   PetscStrallocpy(SNESMonitorFunction,&sctx->funcname);
5007:   /*
5008:      This should work, but it doesn't
5009:   sctx->ctx = ctx;
5010:   mexMakeArrayPersistent(sctx->ctx);
5011:   */
5012:   sctx->ctx = mxDuplicateArray(ctx);
5013:   SNESMonitorSet(snes,SNESMonitor_Matlab,sctx,NULL);
5014:   return(0);
5015: }

5017: #endif