Actual source code: snes.c

petsc-3.11.3 2019-06-26
Report Typos and Errors
  1:  #include <petsc/private/snesimpl.h>
  2:  #include <petscdmshell.h>
  3:  #include <petscdraw.h>
  4:  #include <petscds.h>
  5:  #include <petscdmadaptor.h>
  6:  #include <petscconvest.h>

  8: PetscBool         SNESRegisterAllCalled = PETSC_FALSE;
  9: PetscFunctionList SNESList              = NULL;

 11: /* Logging support */
 12: PetscClassId  SNES_CLASSID, DMSNES_CLASSID;
 13: PetscLogEvent SNES_Solve, SNES_FunctionEval, SNES_JacobianEval, SNES_NGSEval, SNES_NGSFuncEval, SNES_NPCSolve, SNES_ObjectiveEval;

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

 18:    Logically Collective on SNES

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

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

 27:    Level: intermediate

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

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

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

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

 49:    Not Collective

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

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

 57:    Level: intermediate

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

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

 72: /*@
 73:     SNESSetAlwaysComputesFinalResidual - does the SNES always compute the residual at the final solution?

 75:    Logically Collective on SNES

 77:     Input Parameters:
 78: +   snes - the shell SNES
 79: -   flg - is the residual computed?

 81:    Level: advanced

 83: .seealso: SNESGetAlwaysComputesFinalResidual()
 84: @*/
 85: PetscErrorCode  SNESSetAlwaysComputesFinalResidual(SNES snes, PetscBool flg)
 86: {
 89:   snes->alwayscomputesfinalresidual = flg;
 90:   return(0);
 91: }

 93: /*@
 94:     SNESGetAlwaysComputesFinalResidual - does the SNES always compute the residual at the final solution?

 96:    Logically Collective on SNES

 98:     Input Parameter:
 99: .   snes - the shell SNES

101:     Output Parameter:
102: .   flg - is the residual computed?

104:    Level: advanced

106: .seealso: SNESSetAlwaysComputesFinalResidual()
107: @*/
108: PetscErrorCode  SNESGetAlwaysComputesFinalResidual(SNES snes, PetscBool *flg)
109: {
112:   *flg = snes->alwayscomputesfinalresidual;
113:   return(0);
114: }

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

120:    Logically Collective on SNES

122:    Input Parameters:
123: .  snes - the SNES context

125:    Level: advanced

127: .keywords: SNES, view

129: .seealso: SNESCreate(), SNESSetFunction(), SNESFunction
130: @*/
131: PetscErrorCode  SNESSetFunctionDomainError(SNES snes)
132: {
135:   if (snes->errorifnotconverged) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"User code indicates input vector is not in the function domain");
136:   snes->domainerror = PETSC_TRUE;
137:   return(0);
138: }

140: /*@
141:    SNESSetJacobianDomainError - tells SNES that computeJacobian does not make sense any more. For example there is a negative element transformation.

143:    Logically Collective on SNES

145:    Input Parameters:
146: .  snes - the SNES context

148:    Level: advanced

150: .keywords: SNES, view

152: .seealso: SNESCreate(), SNESSetFunction(), SNESFunction(), SNESSetFunctionDomainError()
153: @*/
154: PetscErrorCode SNESSetJacobianDomainError(SNES snes)
155: {
158:   if (snes->errorifnotconverged) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"User code indicates computeJacobian does not make sense");
159:   snes->jacobiandomainerror = PETSC_TRUE;
160:   return(0);
161: }

163: /*@
164:    SNESSetCheckJacobianDomainError - if or not to check jacobian domain error after each Jacobian evaluation. By default, we check Jacobian domain error
165:    in the debug mode, and do not check it in the optimized mode.

167:    Logically Collective on SNES

169:    Input Parameters:
170: .  snes - the SNES context
171: .  flg  - indicates if or not to check jacobian domain error after each Jacobian evaluation

173:    Level: advanced

175: .keywords: SNES, view

177: .seealso: SNESCreate(), SNESSetFunction(), SNESFunction(), SNESSetFunctionDomainError(), SNESGetCheckJacobianDomainError()
178: @*/
179: PetscErrorCode SNESSetCheckJacobianDomainError(SNES snes, PetscBool flg)
180: {
183:   snes->checkjacdomainerror = flg;
184:   return(0);
185: }

187: /*@
188:    SNESGetCheckJacobianDomainError - Get an indicator whether or not we are checking Jacobian domain errors after each Jacobian evaluation.

190:    Logically Collective on SNES

192:    Input Parameters:
193: .  snes - the SNES context

195:    Output Parameters:
196: .  flg  - PETSC_FALSE indicates that we don't check jacobian domain errors after each Jacobian evaluation

198:    Level: advanced

200: .keywords: SNES, view

202: .seealso: SNESCreate(), SNESSetFunction(), SNESFunction(), SNESSetFunctionDomainError(), SNESSetCheckJacobianDomainError()
203: @*/
204: PetscErrorCode SNESGetCheckJacobianDomainError(SNES snes, PetscBool *flg)
205: {
209:   *flg = snes->checkjacdomainerror;
210:   return(0);
211: }

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

216:    Logically Collective on SNES

218:    Input Parameters:
219: .  snes - the SNES context

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

224:    Level: advanced

226: .keywords: SNES, view

228: .seealso: SNESSetFunctionDomainError(), SNESComputeFunction()
229: @*/
230: PetscErrorCode  SNESGetFunctionDomainError(SNES snes, PetscBool *domainerror)
231: {
235:   *domainerror = snes->domainerror;
236:   return(0);
237: }

239: /*@
240:    SNESGetJacobianDomainError - Gets the status of the Jacobian domain error after a call to SNESComputeJacobian;

242:    Logically Collective on SNES

244:    Input Parameters:
245: .  snes - the SNES context

247:    Output Parameters:
248: .  domainerror - Set to PETSC_TRUE if there's a jacobian domain error; PETSC_FALSE otherwise.

250:    Level: advanced

252: .keywords: SNES, view

254: .seealso: SNESSetFunctionDomainError(), SNESComputeFunction(),SNESGetFunctionDomainError()
255: @*/
256: PetscErrorCode SNESGetJacobianDomainError(SNES snes, PetscBool *domainerror)
257: {
261:   *domainerror = snes->jacobiandomainerror;
262:   return(0);
263: }

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

268:   Collective on PetscViewer

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

275:    Level: intermediate

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

280:   Notes for advanced users:
281:   Most users should not need to know the details of the binary storage
282:   format, since SNESLoad() and TSView() completely hide these details.
283:   But for anyone who's interested, the standard binary matrix storage
284:   format is
285: .vb
286:      has not yet been determined
287: .ve

289: .seealso: PetscViewerBinaryOpen(), SNESView(), MatLoad(), VecLoad()
290: @*/
291: PetscErrorCode  SNESLoad(SNES snes, PetscViewer viewer)
292: {
294:   PetscBool      isbinary;
295:   PetscInt       classid;
296:   char           type[256];
297:   KSP            ksp;
298:   DM             dm;
299:   DMSNES         dmsnes;

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

307:   PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
308:   if (classid != SNES_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONG,"Not SNES next in file");
309:   PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
310:   SNESSetType(snes, type);
311:   if (snes->ops->load) {
312:     (*snes->ops->load)(snes,viewer);
313:   }
314:   SNESGetDM(snes,&dm);
315:   DMGetDMSNES(dm,&dmsnes);
316:   DMSNESLoad(dmsnes,viewer);
317:   SNESGetKSP(snes,&ksp);
318:   KSPLoad(ksp,viewer);
319:   return(0);
320: }

322:  #include <petscdraw.h>
323: #if defined(PETSC_HAVE_SAWS)
324:  #include <petscviewersaws.h>
325: #endif

327: /*@C
328:    SNESView - Prints the SNES data structure.

330:    Collective on SNES

332:    Input Parameters:
333: +  SNES - the SNES context
334: -  viewer - visualization context

336:    Options Database Key:
337: .  -snes_view - Calls SNESView() at end of SNESSolve()

339:    Notes:
340:    The available visualization contexts include
341: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
342: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
343:          output where only the first processor opens
344:          the file.  All other processors send their
345:          data to the first processor to print.

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

350:    Level: beginner

352: .keywords: SNES, view

354: .seealso: PetscViewerASCIIOpen()
355: @*/
356: PetscErrorCode  SNESView(SNES snes,PetscViewer viewer)
357: {
358:   SNESKSPEW      *kctx;
360:   KSP            ksp;
361:   SNESLineSearch linesearch;
362:   PetscBool      iascii,isstring,isbinary,isdraw;
363:   DMSNES         dmsnes;
364: #if defined(PETSC_HAVE_SAWS)
365:   PetscBool      issaws;
366: #endif

370:   if (!viewer) {
371:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&viewer);
372:   }

376:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
377:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
378:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
379:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
380: #if defined(PETSC_HAVE_SAWS)
381:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSAWS,&issaws);
382: #endif
383:   if (iascii) {
384:     SNESNormSchedule normschedule;
385:     DM               dm;
386:     PetscErrorCode   (*cJ)(SNES,Vec,Mat,Mat,void*);
387:     void             *ctx;

389:     PetscObjectPrintClassNamePrefixType((PetscObject)snes,viewer);
390:     if (!snes->setupcalled) {
391:       PetscViewerASCIIPrintf(viewer,"  SNES has not been set up so information may be incomplete\n");
392:     }
393:     if (snes->ops->view) {
394:       PetscViewerASCIIPushTab(viewer);
395:       (*snes->ops->view)(snes,viewer);
396:       PetscViewerASCIIPopTab(viewer);
397:     }
398:     PetscViewerASCIIPrintf(viewer,"  maximum iterations=%D, maximum function evaluations=%D\n",snes->max_its,snes->max_funcs);
399:     PetscViewerASCIIPrintf(viewer,"  tolerances: relative=%g, absolute=%g, solution=%g\n",(double)snes->rtol,(double)snes->abstol,(double)snes->stol);
400:     if (snes->usesksp) {
401:       PetscViewerASCIIPrintf(viewer,"  total number of linear solver iterations=%D\n",snes->linear_its);
402:     }
403:     PetscViewerASCIIPrintf(viewer,"  total number of function evaluations=%D\n",snes->nfuncs);
404:     SNESGetNormSchedule(snes, &normschedule);
405:     if (normschedule > 0) {PetscViewerASCIIPrintf(viewer,"  norm schedule %s\n",SNESNormSchedules[normschedule]);}
406:     if (snes->gridsequence) {
407:       PetscViewerASCIIPrintf(viewer,"  total number of grid sequence refinements=%D\n",snes->gridsequence);
408:     }
409:     if (snes->ksp_ewconv) {
410:       kctx = (SNESKSPEW*)snes->kspconvctx;
411:       if (kctx) {
412:         PetscViewerASCIIPrintf(viewer,"  Eisenstat-Walker computation of KSP relative tolerance (version %D)\n",kctx->version);
413:         PetscViewerASCIIPrintf(viewer,"    rtol_0=%g, rtol_max=%g, threshold=%g\n",(double)kctx->rtol_0,(double)kctx->rtol_max,(double)kctx->threshold);
414:         PetscViewerASCIIPrintf(viewer,"    gamma=%g, alpha=%g, alpha2=%g\n",(double)kctx->gamma,(double)kctx->alpha,(double)kctx->alpha2);
415:       }
416:     }
417:     if (snes->lagpreconditioner == -1) {
418:       PetscViewerASCIIPrintf(viewer,"  Preconditioned is never rebuilt\n");
419:     } else if (snes->lagpreconditioner > 1) {
420:       PetscViewerASCIIPrintf(viewer,"  Preconditioned is rebuilt every %D new Jacobians\n",snes->lagpreconditioner);
421:     }
422:     if (snes->lagjacobian == -1) {
423:       PetscViewerASCIIPrintf(viewer,"  Jacobian is never rebuilt\n");
424:     } else if (snes->lagjacobian > 1) {
425:       PetscViewerASCIIPrintf(viewer,"  Jacobian is rebuilt every %D SNES iterations\n",snes->lagjacobian);
426:     }
427:     SNESGetDM(snes,&dm);
428:     DMSNESGetJacobian(dm,&cJ,&ctx);
429:     if (cJ == SNESComputeJacobianDefault) {
430:       PetscViewerASCIIPrintf(viewer,"  Jacobian is built using finite differences one column at a time\n");
431:     } else if (cJ == SNESComputeJacobianDefaultColor) {
432:       PetscViewerASCIIPrintf(viewer,"  Jacobian is built using finite differences with coloring\n");
433:     }
434:   } else if (isstring) {
435:     const char *type;
436:     SNESGetType(snes,&type);
437:     PetscViewerStringSPrintf(viewer," %-3.3s",type);
438:   } else if (isbinary) {
439:     PetscInt    classid = SNES_FILE_CLASSID;
440:     MPI_Comm    comm;
441:     PetscMPIInt rank;
442:     char        type[256];

444:     PetscObjectGetComm((PetscObject)snes,&comm);
445:     MPI_Comm_rank(comm,&rank);
446:     if (!rank) {
447:       PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);
448:       PetscStrncpy(type,((PetscObject)snes)->type_name,sizeof(type));
449:       PetscViewerBinaryWrite(viewer,type,sizeof(type),PETSC_CHAR,PETSC_FALSE);
450:     }
451:     if (snes->ops->view) {
452:       (*snes->ops->view)(snes,viewer);
453:     }
454:   } else if (isdraw) {
455:     PetscDraw draw;
456:     char      str[36];
457:     PetscReal x,y,bottom,h;

459:     PetscViewerDrawGetDraw(viewer,0,&draw);
460:     PetscDrawGetCurrentPoint(draw,&x,&y);
461:     PetscStrncpy(str,"SNES: ",sizeof(str));
462:     PetscStrlcat(str,((PetscObject)snes)->type_name,sizeof(str));
463:     PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_BLUE,PETSC_DRAW_BLACK,str,NULL,&h);
464:     bottom = y - h;
465:     PetscDrawPushCurrentPoint(draw,x,bottom);
466:     if (snes->ops->view) {
467:       (*snes->ops->view)(snes,viewer);
468:     }
469: #if defined(PETSC_HAVE_SAWS)
470:   } else if (issaws) {
471:     PetscMPIInt rank;
472:     const char *name;

474:     PetscObjectGetName((PetscObject)snes,&name);
475:     MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
476:     if (!((PetscObject)snes)->amsmem && !rank) {
477:       char       dir[1024];

479:       PetscObjectViewSAWs((PetscObject)snes,viewer);
480:       PetscSNPrintf(dir,1024,"/PETSc/Objects/%s/its",name);
481:       PetscStackCallSAWs(SAWs_Register,(dir,&snes->iter,1,SAWs_READ,SAWs_INT));
482:       if (!snes->conv_hist) {
483:         SNESSetConvergenceHistory(snes,NULL,NULL,PETSC_DECIDE,PETSC_TRUE);
484:       }
485:       PetscSNPrintf(dir,1024,"/PETSc/Objects/%s/conv_hist",name);
486:       PetscStackCallSAWs(SAWs_Register,(dir,snes->conv_hist,10,SAWs_READ,SAWs_DOUBLE));
487:     }
488: #endif
489:   }
490:   if (snes->linesearch) {
491:     SNESGetLineSearch(snes, &linesearch);
492:     PetscViewerASCIIPushTab(viewer);
493:     SNESLineSearchView(linesearch, viewer);
494:     PetscViewerASCIIPopTab(viewer);
495:   }
496:   if (snes->npc && snes->usesnpc) {
497:     PetscViewerASCIIPushTab(viewer);
498:     SNESView(snes->npc, viewer);
499:     PetscViewerASCIIPopTab(viewer);
500:   }
501:   PetscViewerASCIIPushTab(viewer);
502:   DMGetDMSNES(snes->dm,&dmsnes);
503:   DMSNESView(dmsnes, viewer);
504:   PetscViewerASCIIPopTab(viewer);
505:   if (snes->usesksp) {
506:     SNESGetKSP(snes,&ksp);
507:     PetscViewerASCIIPushTab(viewer);
508:     KSPView(ksp,viewer);
509:     PetscViewerASCIIPopTab(viewer);
510:   }
511:   if (isdraw) {
512:     PetscDraw draw;
513:     PetscViewerDrawGetDraw(viewer,0,&draw);
514:     PetscDrawPopCurrentPoint(draw);
515:   }
516:   return(0);
517: }

519: /*
520:   We retain a list of functions that also take SNES command
521:   line options. These are called at the end SNESSetFromOptions()
522: */
523: #define MAXSETFROMOPTIONS 5
524: static PetscInt numberofsetfromoptions;
525: static PetscErrorCode (*othersetfromoptions[MAXSETFROMOPTIONS])(SNES);

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

530:   Not Collective

532:   Input Parameter:
533: . snescheck - function that checks for options

535:   Level: developer

537: .seealso: SNESSetFromOptions()
538: @*/
539: PetscErrorCode  SNESAddOptionsChecker(PetscErrorCode (*snescheck)(SNES))
540: {
542:   if (numberofsetfromoptions >= MAXSETFROMOPTIONS) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Too many options checkers, only %D allowed", MAXSETFROMOPTIONS);
543:   othersetfromoptions[numberofsetfromoptions++] = snescheck;
544:   return(0);
545: }

547: PETSC_INTERN PetscErrorCode SNESDefaultMatrixFreeCreate2(SNES,Vec,Mat*);

549: static PetscErrorCode SNESSetUpMatrixFree_Private(SNES snes, PetscBool hasOperator, PetscInt version)
550: {
551:   Mat            J;
552:   KSP            ksp;
553:   PC             pc;
554:   PetscBool      match;
556:   MatNullSpace   nullsp;


561:   if (!snes->vec_func && (snes->jacobian || snes->jacobian_pre)) {
562:     Mat A = snes->jacobian, B = snes->jacobian_pre;
563:     MatCreateVecs(A ? A : B, NULL,&snes->vec_func);
564:   }

566:   if (version == 1) {
567:     MatCreateSNESMF(snes,&J);
568:     MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);
569:     MatSetFromOptions(J);
570:   } else if (version == 2) {
571:     if (!snes->vec_func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"SNESSetFunction() must be called first");
572: #if !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128) && !defined(PETSC_USE_REAL___FP16)
573:     SNESDefaultMatrixFreeCreate2(snes,snes->vec_func,&J);
574: #else
575:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "matrix-free operator rutines (version 2)");
576: #endif
577:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "matrix-free operator rutines, only version 1 and 2");

579:   /* attach any user provided null space that was on Amat to the newly created matrix free matrix */
580:   if (snes->jacobian) {
581:     MatGetNullSpace(snes->jacobian,&nullsp);
582:     if (nullsp) {
583:       MatSetNullSpace(J,nullsp);
584:     }
585:   }

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

590:     /* This version replaces the user provided Jacobian matrix with a
591:        matrix-free version but still employs the user-provided preconditioner matrix. */
592:     SNESSetJacobian(snes,J,0,0,0);
593:   } else {
594:     /* This version replaces both the user-provided Jacobian and the user-
595:      provided preconditioner Jacobian with the default matrix free version. */
596:     if ((snes->npcside== PC_LEFT) && snes->npc) {
597:       if (!snes->jacobian){SNESSetJacobian(snes,J,0,0,0);}
598:     } else {
599:       SNESSetJacobian(snes,J,J,MatMFFDComputeJacobian,0);
600:     }
601:     /* Force no preconditioner */
602:     SNESGetKSP(snes,&ksp);
603:     KSPGetPC(ksp,&pc);
604:     PetscObjectTypeCompare((PetscObject)pc,PCSHELL,&match);
605:     if (!match) {
606:       PetscInfo(snes,"Setting default matrix-free preconditioner routines\nThat is no preconditioner is being used\n");
607:       PCSetType(pc,PCNONE);
608:     }
609:   }
610:   MatDestroy(&J);
611:   return(0);
612: }

614: static PetscErrorCode DMRestrictHook_SNESVecSol(DM dmfine,Mat Restrict,Vec Rscale,Mat Inject,DM dmcoarse,void *ctx)
615: {
616:   SNES           snes = (SNES)ctx;
618:   Vec            Xfine,Xfine_named = NULL,Xcoarse;

621:   if (PetscLogPrintInfo) {
622:     PetscInt finelevel,coarselevel,fineclevel,coarseclevel;
623:     DMGetRefineLevel(dmfine,&finelevel);
624:     DMGetCoarsenLevel(dmfine,&fineclevel);
625:     DMGetRefineLevel(dmcoarse,&coarselevel);
626:     DMGetCoarsenLevel(dmcoarse,&coarseclevel);
627:     PetscInfo4(dmfine,"Restricting SNES solution vector from level %D-%D to level %D-%D\n",finelevel,fineclevel,coarselevel,coarseclevel);
628:   }
629:   if (dmfine == snes->dm) Xfine = snes->vec_sol;
630:   else {
631:     DMGetNamedGlobalVector(dmfine,"SNESVecSol",&Xfine_named);
632:     Xfine = Xfine_named;
633:   }
634:   DMGetNamedGlobalVector(dmcoarse,"SNESVecSol",&Xcoarse);
635:   if (Inject) {
636:     MatRestrict(Inject,Xfine,Xcoarse);
637:   } else {
638:     MatRestrict(Restrict,Xfine,Xcoarse);
639:     VecPointwiseMult(Xcoarse,Xcoarse,Rscale);
640:   }
641:   DMRestoreNamedGlobalVector(dmcoarse,"SNESVecSol",&Xcoarse);
642:   if (Xfine_named) {DMRestoreNamedGlobalVector(dmfine,"SNESVecSol",&Xfine_named);}
643:   return(0);
644: }

646: static PetscErrorCode DMCoarsenHook_SNESVecSol(DM dm,DM dmc,void *ctx)
647: {

651:   DMCoarsenHookAdd(dmc,DMCoarsenHook_SNESVecSol,DMRestrictHook_SNESVecSol,ctx);
652:   return(0);
653: }

655: /* This may be called to rediscretize the operator on levels of linear multigrid. The DM shuffle is so the user can
656:  * safely call SNESGetDM() in their residual evaluation routine. */
657: static PetscErrorCode KSPComputeOperators_SNES(KSP ksp,Mat A,Mat B,void *ctx)
658: {
659:   SNES           snes = (SNES)ctx;
661:   Vec            X,Xnamed = NULL;
662:   DM             dmsave;
663:   void           *ctxsave;
664:   PetscErrorCode (*jac)(SNES,Vec,Mat,Mat,void*) = NULL;

667:   dmsave = snes->dm;
668:   KSPGetDM(ksp,&snes->dm);
669:   if (dmsave == snes->dm) X = snes->vec_sol; /* We are on the finest level */
670:   else {                                     /* We are on a coarser level, this vec was initialized using a DM restrict hook */
671:     DMGetNamedGlobalVector(snes->dm,"SNESVecSol",&Xnamed);
672:     X    = Xnamed;
673:     SNESGetJacobian(snes,NULL,NULL,&jac,&ctxsave);
674:     /* If the DM's don't match up, the MatFDColoring context needed for the jacobian won't match up either -- fixit. */
675:     if (jac == SNESComputeJacobianDefaultColor) {
676:       SNESSetJacobian(snes,NULL,NULL,SNESComputeJacobianDefaultColor,0);
677:     }
678:   }
679:   /* Make sure KSP DM has the Jacobian computation routine */
680:   {
681:     DMSNES sdm;

683:     DMGetDMSNES(snes->dm, &sdm);
684:     if (!sdm->ops->computejacobian) {
685:       DMCopyDMSNES(dmsave, snes->dm);
686:     }
687:   }
688:   /* Compute the operators */
689:   SNESComputeJacobian(snes,X,A,B);
690:   /* Put the previous context back */
691:   if (snes->dm != dmsave && jac == SNESComputeJacobianDefaultColor) {
692:     SNESSetJacobian(snes,NULL,NULL,jac,ctxsave);
693:   }

695:   if (Xnamed) {DMRestoreNamedGlobalVector(snes->dm,"SNESVecSol",&Xnamed);}
696:   snes->dm = dmsave;
697:   return(0);
698: }

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

703:    Collective

705:    Input Arguments:
706: .  snes - snes to configure

708:    Level: developer

710: .seealso: SNESSetUp()
711: @*/
712: PetscErrorCode SNESSetUpMatrices(SNES snes)
713: {
715:   DM             dm;
716:   DMSNES         sdm;

719:   SNESGetDM(snes,&dm);
720:   DMGetDMSNES(dm,&sdm);
721:   if (!sdm->ops->computejacobian) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"DMSNES not properly configured");
722:   else if (!snes->jacobian && snes->mf) {
723:     Mat  J;
724:     void *functx;
725:     MatCreateSNESMF(snes,&J);
726:     MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);
727:     MatSetFromOptions(J);
728:     SNESGetFunction(snes,NULL,NULL,&functx);
729:     SNESSetJacobian(snes,J,J,0,0);
730:     MatDestroy(&J);
731:   } else if (snes->mf_operator && !snes->jacobian_pre && !snes->jacobian) {
732:     Mat J,B;
733:     MatCreateSNESMF(snes,&J);
734:     MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);
735:     MatSetFromOptions(J);
736:     DMCreateMatrix(snes->dm,&B);
737:     /* sdm->computejacobian was already set to reach here */
738:     SNESSetJacobian(snes,J,B,NULL,NULL);
739:     MatDestroy(&J);
740:     MatDestroy(&B);
741:   } else if (!snes->jacobian_pre) {
742:     PetscErrorCode (*nspconstr)(DM, PetscInt, MatNullSpace *);
743:     PetscDS          prob;
744:     Mat              J, B;
745:     MatNullSpace     nullspace = NULL;
746:     PetscBool        hasPrec   = PETSC_FALSE;
747:     PetscInt         Nf;

749:     J    = snes->jacobian;
750:     DMGetDS(dm, &prob);
751:     if (prob) {PetscDSHasJacobianPreconditioner(prob, &hasPrec);}
752:     if (J)            {PetscObjectReference((PetscObject) J);}
753:     else if (hasPrec) {DMCreateMatrix(snes->dm, &J);}
754:     DMCreateMatrix(snes->dm, &B);
755:     PetscDSGetNumFields(prob, &Nf);
756:     DMGetNullSpaceConstructor(snes->dm, Nf, &nspconstr);
757:     if (nspconstr) (*nspconstr)(snes->dm, -1, &nullspace);
758:     MatSetNullSpace(B, nullspace);
759:     MatNullSpaceDestroy(&nullspace);
760:     SNESSetJacobian(snes, J ? J : B, B, NULL, NULL);
761:     MatDestroy(&J);
762:     MatDestroy(&B);
763:   }
764:   {
765:     KSP ksp;
766:     SNESGetKSP(snes,&ksp);
767:     KSPSetComputeOperators(ksp,KSPComputeOperators_SNES,snes);
768:     DMCoarsenHookAdd(snes->dm,DMCoarsenHook_SNESVecSol,DMRestrictHook_SNESVecSol,snes);
769:   }
770:   return(0);
771: }

773: /*@C
774:    SNESMonitorSetFromOptions - Sets a monitor function and viewer appropriate for the type indicated by the user

776:    Collective on SNES

778:    Input Parameters:
779: +  snes - SNES object you wish to monitor
780: .  name - the monitor type one is seeking
781: .  help - message indicating what monitoring is done
782: .  manual - manual page for the monitor
783: .  monitor - the monitor function
784: -  monitorsetup - a function that is called once ONLY if the user selected this monitor that may set additional features of the SNES or PetscViewer objects

786:    Level: developer

788: .seealso: PetscOptionsGetViewer(), PetscOptionsGetReal(), PetscOptionsHasName(), PetscOptionsGetString(),
789:           PetscOptionsGetIntArray(), PetscOptionsGetRealArray(), PetscOptionsBool()
790:           PetscOptionsInt(), PetscOptionsString(), PetscOptionsReal(), PetscOptionsBool(),
791:           PetscOptionsName(), PetscOptionsBegin(), PetscOptionsEnd(), PetscOptionsHead(),
792:           PetscOptionsStringArray(),PetscOptionsRealArray(), PetscOptionsScalar(),
793:           PetscOptionsBoolGroupBegin(), PetscOptionsBoolGroup(), PetscOptionsBoolGroupEnd(),
794:           PetscOptionsFList(), PetscOptionsEList()
795: @*/
796: PetscErrorCode  SNESMonitorSetFromOptions(SNES snes,const char name[],const char help[], const char manual[],PetscErrorCode (*monitor)(SNES,PetscInt,PetscReal,PetscViewerAndFormat*),PetscErrorCode (*monitorsetup)(SNES,PetscViewerAndFormat*))
797: {
798:   PetscErrorCode    ierr;
799:   PetscViewer       viewer;
800:   PetscViewerFormat format;
801:   PetscBool         flg;

804:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject) snes)->options,((PetscObject)snes)->prefix,name,&viewer,&format,&flg);
805:   if (flg) {
806:     PetscViewerAndFormat *vf;
807:     PetscViewerAndFormatCreate(viewer,format,&vf);
808:     PetscObjectDereference((PetscObject)viewer);
809:     if (monitorsetup) {
810:       (*monitorsetup)(snes,vf);
811:     }
812:     SNESMonitorSet(snes,(PetscErrorCode (*)(SNES,PetscInt,PetscReal,void*))monitor,vf,(PetscErrorCode (*)(void**))PetscViewerAndFormatDestroy);
813:   }
814:   return(0);
815: }

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

820:    Collective on SNES

822:    Input Parameter:
823: .  snes - the SNES context

825:    Options Database Keys:
826: +  -snes_type <type> - newtonls, newtontr, ngmres, ncg, nrichardson, qn, vi, fas, SNESType for complete list
827: .  -snes_stol - convergence tolerance in terms of the norm
828:                 of the change in the solution between steps
829: .  -snes_atol <abstol> - absolute tolerance of residual norm
830: .  -snes_rtol <rtol> - relative decrease in tolerance norm from initial
831: .  -snes_divergence_tolerance <divtol> - if the residual goes above divtol*rnorm0, exit with divergence
832: .  -snes_force_iteration <force> - force SNESSolve() to take at least one iteration
833: .  -snes_max_it <max_it> - maximum number of iterations
834: .  -snes_max_funcs <max_funcs> - maximum number of function evaluations
835: .  -snes_max_fail <max_fail> - maximum number of line search failures allowed before stopping, default is none
836: .  -snes_max_linear_solve_fail - number of linear solver failures before SNESSolve() stops
837: .  -snes_lag_preconditioner <lag> - how often preconditioner is rebuilt (use -1 to never rebuild)
838: .  -snes_lag_jacobian <lag> - how often Jacobian is rebuilt (use -1 to never rebuild)
839: .  -snes_trtol <trtol> - trust region tolerance
840: .  -snes_no_convergence_test - skip convergence test in nonlinear
841:                                solver; hence iterations will continue until max_it
842:                                or some other criterion is reached. Saves expense
843:                                of convergence test
844: .  -snes_monitor [ascii][:filename][:viewer format] - prints residual norm at each iteration. if no filename given prints to stdout
845: .  -snes_monitor_solution [ascii binary draw][:filename][:viewer format] - plots solution at each iteration
846: .  -snes_monitor_residual [ascii binary draw][:filename][:viewer format] - plots residual (not its norm) at each iteration
847: .  -snes_monitor_solution_update [ascii binary draw][:filename][:viewer format] - plots update to solution at each iteration
848: .  -snes_monitor_lg_residualnorm - plots residual norm at each iteration
849: .  -snes_monitor_lg_range - plots residual norm at each iteration
850: .  -snes_fd - use finite differences to compute Jacobian; very slow, only for testing
851: .  -snes_fd_color - use finite differences with coloring to compute Jacobian
852: .  -snes_mf_ksp_monitor - if using matrix-free multiply then print h at each KSP iteration
853: -  -snes_converged_reason - print the reason for convergence/divergence after each solve

855:     Options Database for Eisenstat-Walker method:
856: +  -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence
857: .  -snes_ksp_ew_version ver - version of  Eisenstat-Walker method
858: .  -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0
859: .  -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax
860: .  -snes_ksp_ew_gamma <gamma> - Sets gamma
861: .  -snes_ksp_ew_alpha <alpha> - Sets alpha
862: .  -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2
863: -  -snes_ksp_ew_threshold <threshold> - Sets threshold

865:    Notes:
866:    To see all options, run your program with the -help option or consult
867:    Users-Manual: ch_snes

869:    Level: beginner

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

873: .seealso: SNESSetOptionsPrefix(), SNESResetFromOptions()
874: @*/
875: PetscErrorCode  SNESSetFromOptions(SNES snes)
876: {
877:   PetscBool      flg,pcset,persist,set;
878:   PetscInt       i,indx,lag,grids;
879:   const char     *deft        = SNESNEWTONLS;
880:   const char     *convtests[] = {"default","skip"};
881:   SNESKSPEW      *kctx        = NULL;
882:   char           type[256], monfilename[PETSC_MAX_PATH_LEN];
884:   PCSide         pcside;
885:   const char     *optionsprefix;

889:   SNESRegisterAll();
890:   PetscObjectOptionsBegin((PetscObject)snes);
891:   if (((PetscObject)snes)->type_name) deft = ((PetscObject)snes)->type_name;
892:   PetscOptionsFList("-snes_type","Nonlinear solver method","SNESSetType",SNESList,deft,type,256,&flg);
893:   if (flg) {
894:     SNESSetType(snes,type);
895:   } else if (!((PetscObject)snes)->type_name) {
896:     SNESSetType(snes,deft);
897:   }
898:   PetscOptionsReal("-snes_stol","Stop if step length less than","SNESSetTolerances",snes->stol,&snes->stol,NULL);
899:   PetscOptionsReal("-snes_atol","Stop if function norm less than","SNESSetTolerances",snes->abstol,&snes->abstol,NULL);

901:   PetscOptionsReal("-snes_rtol","Stop if decrease in function norm less than","SNESSetTolerances",snes->rtol,&snes->rtol,NULL);
902:   PetscOptionsReal("-snes_divergence_tolerance","Stop if residual norm increases by this factor","SNESSetDivergenceTolerance",snes->divtol,&snes->divtol,NULL);
903:   PetscOptionsInt("-snes_max_it","Maximum iterations","SNESSetTolerances",snes->max_its,&snes->max_its,NULL);
904:   PetscOptionsInt("-snes_max_funcs","Maximum function evaluations","SNESSetTolerances",snes->max_funcs,&snes->max_funcs,NULL);
905:   PetscOptionsInt("-snes_max_fail","Maximum nonlinear step failures","SNESSetMaxNonlinearStepFailures",snes->maxFailures,&snes->maxFailures,NULL);
906:   PetscOptionsInt("-snes_max_linear_solve_fail","Maximum failures in linear solves allowed","SNESSetMaxLinearSolveFailures",snes->maxLinearSolveFailures,&snes->maxLinearSolveFailures,NULL);
907:   PetscOptionsBool("-snes_error_if_not_converged","Generate error if solver does not converge","SNESSetErrorIfNotConverged",snes->errorifnotconverged,&snes->errorifnotconverged,NULL);
908:   PetscOptionsBool("-snes_force_iteration","Force SNESSolve() to take at least one iteration","SNESSetForceIteration",snes->forceiteration,&snes->forceiteration,NULL);
909:   PetscOptionsBool("-snes_check_jacobian_domain_error","Check Jacobian domain error after Jacobian evaluation","SNESCheckJacobianDomainError",snes->checkjacdomainerror,&snes->checkjacdomainerror,NULL);

911:   PetscOptionsInt("-snes_lag_preconditioner","How often to rebuild preconditioner","SNESSetLagPreconditioner",snes->lagpreconditioner,&lag,&flg);
912:   if (flg) {
913:     SNESSetLagPreconditioner(snes,lag);
914:   }
915:   PetscOptionsBool("-snes_lag_preconditioner_persists","Preconditioner lagging through multiple solves","SNESSetLagPreconditionerPersists",snes->lagjac_persist,&persist,&flg);
916:   if (flg) {
917:     SNESSetLagPreconditionerPersists(snes,persist);
918:   }
919:   PetscOptionsInt("-snes_lag_jacobian","How often to rebuild Jacobian","SNESSetLagJacobian",snes->lagjacobian,&lag,&flg);
920:   if (flg) {
921:     SNESSetLagJacobian(snes,lag);
922:   }
923:   PetscOptionsBool("-snes_lag_jacobian_persists","Jacobian lagging through multiple solves","SNESSetLagJacobianPersists",snes->lagjac_persist,&persist,&flg);
924:   if (flg) {
925:     SNESSetLagJacobianPersists(snes,persist);
926:   }

928:   PetscOptionsInt("-snes_grid_sequence","Use grid sequencing to generate initial guess","SNESSetGridSequence",snes->gridsequence,&grids,&flg);
929:   if (flg) {
930:     SNESSetGridSequence(snes,grids);
931:   }

933:   PetscOptionsEList("-snes_convergence_test","Convergence test","SNESSetConvergenceTest",convtests,2,"default",&indx,&flg);
934:   if (flg) {
935:     switch (indx) {
936:     case 0: SNESSetConvergenceTest(snes,SNESConvergedDefault,NULL,NULL); break;
937:     case 1: SNESSetConvergenceTest(snes,SNESConvergedSkip,NULL,NULL);    break;
938:     }
939:   }

941:   PetscOptionsEList("-snes_norm_schedule","SNES Norm schedule","SNESSetNormSchedule",SNESNormSchedules,5,"function",&indx,&flg);
942:   if (flg) { SNESSetNormSchedule(snes,(SNESNormSchedule)indx); }

944:   PetscOptionsEList("-snes_function_type","SNES Norm schedule","SNESSetFunctionType",SNESFunctionTypes,2,"unpreconditioned",&indx,&flg);
945:   if (flg) { SNESSetFunctionType(snes,(SNESFunctionType)indx); }

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

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

951:   PetscOptionsInt("-snes_ksp_ew_version","Version 1, 2 or 3","SNESKSPSetParametersEW",kctx->version,&kctx->version,NULL);
952:   PetscOptionsReal("-snes_ksp_ew_rtol0","0 <= rtol0 < 1","SNESKSPSetParametersEW",kctx->rtol_0,&kctx->rtol_0,NULL);
953:   PetscOptionsReal("-snes_ksp_ew_rtolmax","0 <= rtolmax < 1","SNESKSPSetParametersEW",kctx->rtol_max,&kctx->rtol_max,NULL);
954:   PetscOptionsReal("-snes_ksp_ew_gamma","0 <= gamma <= 1","SNESKSPSetParametersEW",kctx->gamma,&kctx->gamma,NULL);
955:   PetscOptionsReal("-snes_ksp_ew_alpha","1 < alpha <= 2","SNESKSPSetParametersEW",kctx->alpha,&kctx->alpha,NULL);
956:   PetscOptionsReal("-snes_ksp_ew_alpha2","alpha2","SNESKSPSetParametersEW",kctx->alpha2,&kctx->alpha2,NULL);
957:   PetscOptionsReal("-snes_ksp_ew_threshold","0 < threshold < 1","SNESKSPSetParametersEW",kctx->threshold,&kctx->threshold,NULL);

959:   flg  = PETSC_FALSE;
960:   PetscOptionsBool("-snes_monitor_cancel","Remove all monitors","SNESMonitorCancel",flg,&flg,&set);
961:   if (set && flg) {SNESMonitorCancel(snes);}

963:   SNESMonitorSetFromOptions(snes,"-snes_monitor","Monitor norm of function","SNESMonitorDefault",SNESMonitorDefault,NULL);
964:   SNESMonitorSetFromOptions(snes,"-snes_monitor_short","Monitor norm of function with fewer digits","SNESMonitorDefaultShort",SNESMonitorDefaultShort,NULL);
965:   SNESMonitorSetFromOptions(snes,"-snes_monitor_range","Monitor range of elements of function","SNESMonitorRange",SNESMonitorRange,NULL);

967:   SNESMonitorSetFromOptions(snes,"-snes_monitor_ratio","Monitor ratios of the norm of function for consecutive steps","SNESMonitorRatio",SNESMonitorRatio,SNESMonitorRatioSetUp);
968:   SNESMonitorSetFromOptions(snes,"-snes_monitor_field","Monitor norm of function (split into fields)","SNESMonitorDefaultField",SNESMonitorDefaultField,NULL);
969:   SNESMonitorSetFromOptions(snes,"-snes_monitor_solution","View solution at each iteration","SNESMonitorSolution",SNESMonitorSolution,NULL);
970:   SNESMonitorSetFromOptions(snes,"-snes_monitor_solution_update","View correction at each iteration","SNESMonitorSolutionUpdate",SNESMonitorSolutionUpdate,NULL);
971:   SNESMonitorSetFromOptions(snes,"-snes_monitor_residual","View residual at each iteration","SNESMonitorResidual",SNESMonitorResidual,NULL);
972:   SNESMonitorSetFromOptions(snes,"-snes_monitor_jacupdate_spectrum","Print the change in the spectrum of the Jacobian","SNESMonitorJacUpdateSpectrum",SNESMonitorJacUpdateSpectrum,NULL);
973:   SNESMonitorSetFromOptions(snes,"-snes_monitor_fields","Monitor norm of function per field","SNESMonitorSet",SNESMonitorFields,NULL);

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


979:   flg  = PETSC_FALSE;
980:   PetscOptionsBool("-snes_monitor_lg_residualnorm","Plot function norm at each iteration","SNESMonitorLGResidualNorm",flg,&flg,NULL);
981:   if (flg) {
982:     PetscDrawLG ctx;

984:     SNESMonitorLGCreate(PetscObjectComm((PetscObject)snes),NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,400,300,&ctx);
985:     SNESMonitorSet(snes,SNESMonitorLGResidualNorm,ctx,(PetscErrorCode (*)(void**))PetscDrawLGDestroy);
986:   }
987:   flg  = PETSC_FALSE;
988:   PetscOptionsBool("-snes_monitor_lg_range","Plot function range at each iteration","SNESMonitorLGRange",flg,&flg,NULL);
989:   if (flg) {
990:     PetscViewer ctx;

992:     PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes),NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,400,300,&ctx);
993:     SNESMonitorSet(snes,SNESMonitorLGRange,ctx,(PetscErrorCode (*)(void**))PetscViewerDestroy);
994:   }



998:   flg  = PETSC_FALSE;
999:   PetscOptionsBool("-snes_fd","Use finite differences (slow) to compute Jacobian","SNESComputeJacobianDefault",flg,&flg,NULL);
1000:   if (flg) {
1001:     void    *functx;
1002:     DM      dm;
1003:     DMSNES  sdm;
1004:     SNESGetDM(snes,&dm);
1005:     DMGetDMSNES(dm,&sdm);
1006:     sdm->jacobianctx = NULL;
1007:     SNESGetFunction(snes,NULL,NULL,&functx);
1008:     SNESSetJacobian(snes,snes->jacobian,snes->jacobian_pre,SNESComputeJacobianDefault,functx);
1009:     PetscInfo(snes,"Setting default finite difference Jacobian matrix\n");
1010:   }

1012:   flg  = PETSC_FALSE;
1013:   PetscOptionsBool("-snes_fd_function","Use finite differences (slow) to compute function from user objective","SNESObjectiveComputeFunctionDefaultFD",flg,&flg,NULL);
1014:   if (flg) {
1015:     SNESSetFunction(snes,NULL,SNESObjectiveComputeFunctionDefaultFD,NULL);
1016:   }

1018:   flg  = PETSC_FALSE;
1019:   PetscOptionsBool("-snes_fd_color","Use finite differences with coloring to compute Jacobian","SNESComputeJacobianDefaultColor",flg,&flg,NULL);
1020:   if (flg) {
1021:     DM             dm;
1022:     DMSNES         sdm;
1023:     SNESGetDM(snes,&dm);
1024:     DMGetDMSNES(dm,&sdm);
1025:     sdm->jacobianctx = NULL;
1026:     SNESSetJacobian(snes,snes->jacobian,snes->jacobian_pre,SNESComputeJacobianDefaultColor,0);
1027:     PetscInfo(snes,"Setting default finite difference coloring Jacobian matrix\n");
1028:   }

1030:   flg  = PETSC_FALSE;
1031:   PetscOptionsBool("-snes_mf_operator","Use a Matrix-Free Jacobian with user-provided preconditioner matrix","SNESSetUseMatrixFree",PETSC_FALSE,&snes->mf_operator,&flg);
1032:   if (flg && snes->mf_operator) {
1033:     snes->mf_operator = PETSC_TRUE;
1034:     snes->mf          = PETSC_TRUE;
1035:   }
1036:   flg  = PETSC_FALSE;
1037:   PetscOptionsBool("-snes_mf","Use a Matrix-Free Jacobian with no preconditioner matrix","SNESSetUseMatrixFree",PETSC_FALSE,&snes->mf,&flg);
1038:   if (!flg && snes->mf_operator) snes->mf = PETSC_TRUE;
1039:   PetscOptionsInt("-snes_mf_version","Matrix-Free routines version 1 or 2","None",snes->mf_version,&snes->mf_version,0);

1041:   flg  = PETSC_FALSE;
1042:   SNESGetNPCSide(snes,&pcside);
1043:   PetscOptionsEnum("-snes_npc_side","SNES nonlinear preconditioner side","SNESSetNPCSide",PCSides,(PetscEnum)pcside,(PetscEnum*)&pcside,&flg);
1044:   if (flg) {SNESSetNPCSide(snes,pcside);}

1046: #if defined(PETSC_HAVE_SAWS)
1047:   /*
1048:     Publish convergence information using SAWs
1049:   */
1050:   flg  = PETSC_FALSE;
1051:   PetscOptionsBool("-snes_monitor_saws","Publish SNES progress using SAWs","SNESMonitorSet",flg,&flg,NULL);
1052:   if (flg) {
1053:     void *ctx;
1054:     SNESMonitorSAWsCreate(snes,&ctx);
1055:     SNESMonitorSet(snes,SNESMonitorSAWs,ctx,SNESMonitorSAWsDestroy);
1056:   }
1057: #endif
1058: #if defined(PETSC_HAVE_SAWS)
1059:   {
1060:   PetscBool set;
1061:   flg  = PETSC_FALSE;
1062:   PetscOptionsBool("-snes_saws_block","Block for SAWs at end of SNESSolve","PetscObjectSAWsBlock",((PetscObject)snes)->amspublishblock,&flg,&set);
1063:   if (set) {
1064:     PetscObjectSAWsSetBlock((PetscObject)snes,flg);
1065:   }
1066:   }
1067: #endif

1069:   for (i = 0; i < numberofsetfromoptions; i++) {
1070:     (*othersetfromoptions[i])(snes);
1071:   }

1073:   if (snes->ops->setfromoptions) {
1074:     (*snes->ops->setfromoptions)(PetscOptionsObject,snes);
1075:   }

1077:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
1078:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)snes);
1079:   PetscOptionsEnd();

1081:   if (!snes->linesearch) {
1082:     SNESGetLineSearch(snes, &snes->linesearch);
1083:   }
1084:   SNESLineSearchSetFromOptions(snes->linesearch);

1086:   if (!snes->ksp) {SNESGetKSP(snes,&snes->ksp);}
1087:   KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);
1088:   KSPSetFromOptions(snes->ksp);

1090:   /* if someone has set the SNES NPC type, create it. */
1091:   SNESGetOptionsPrefix(snes, &optionsprefix);
1092:   PetscOptionsHasName(((PetscObject)snes)->options,optionsprefix, "-npc_snes_type", &pcset);
1093:   if (pcset && (!snes->npc)) {
1094:     SNESGetNPC(snes, &snes->npc);
1095:   }
1096:   snes->setfromoptionscalled++;
1097:   return(0);
1098: }

1100: /*@
1101:    SNESResetFromOptions - Sets various SNES and KSP parameters from user options ONLY if the SNES was previously set from options

1103:    Collective on SNES

1105:    Input Parameter:
1106: .  snes - the SNES context

1108:    Level: beginner

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

1112: .seealso: SNESSetFromOptions(), SNESSetOptionsPrefix()
1113: @*/
1114: PetscErrorCode SNESResetFromOptions(SNES snes)
1115: {

1119:   if (snes->setfromoptionscalled) {SNESSetFromOptions(snes);}
1120:   return(0);
1121: }

1123: /*@C
1124:    SNESSetComputeApplicationContext - Sets an optional function to compute a user-defined context for
1125:    the nonlinear solvers.

1127:    Logically Collective on SNES

1129:    Input Parameters:
1130: +  snes - the SNES context
1131: .  compute - function to compute the context
1132: -  destroy - function to destroy the context

1134:    Level: intermediate

1136:    Notes:
1137:    This function is currently not available from Fortran.

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

1141: .seealso: SNESGetApplicationContext(), SNESSetComputeApplicationContext(), SNESGetApplicationContext()
1142: @*/
1143: PetscErrorCode  SNESSetComputeApplicationContext(SNES snes,PetscErrorCode (*compute)(SNES,void**),PetscErrorCode (*destroy)(void**))
1144: {
1147:   snes->ops->usercompute = compute;
1148:   snes->ops->userdestroy = destroy;
1149:   return(0);
1150: }

1152: /*@
1153:    SNESSetApplicationContext - Sets the optional user-defined context for
1154:    the nonlinear solvers.

1156:    Logically Collective on SNES

1158:    Input Parameters:
1159: +  snes - the SNES context
1160: -  usrP - optional user context

1162:    Level: intermediate

1164:    Fortran Notes:
1165:     To use this from Fortran you must write a Fortran interface definition for this
1166:     function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.

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

1170: .seealso: SNESGetApplicationContext()
1171: @*/
1172: PetscErrorCode  SNESSetApplicationContext(SNES snes,void *usrP)
1173: {
1175:   KSP            ksp;

1179:   SNESGetKSP(snes,&ksp);
1180:   KSPSetApplicationContext(ksp,usrP);
1181:   snes->user = usrP;
1182:   return(0);
1183: }

1185: /*@
1186:    SNESGetApplicationContext - Gets the user-defined context for the
1187:    nonlinear solvers.

1189:    Not Collective

1191:    Input Parameter:
1192: .  snes - SNES context

1194:    Output Parameter:
1195: .  usrP - user context

1197:    Fortran Notes:
1198:     To use this from Fortran you must write a Fortran interface definition for this
1199:     function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.

1201:    Level: intermediate

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

1205: .seealso: SNESSetApplicationContext()
1206: @*/
1207: PetscErrorCode  SNESGetApplicationContext(SNES snes,void *usrP)
1208: {
1211:   *(void**)usrP = snes->user;
1212:   return(0);
1213: }

1215: /*@
1216:    SNESSetUseMatrixFree - indicates that SNES should use matrix free finite difference matrix vector products internally to apply
1217:                           the Jacobian.

1219:    Collective on SNES

1221:    Input Parameters:
1222: +  snes - SNES context
1223: .  mf - use matrix-free for both the Amat and Pmat used by SNESSetJacobian(), both the Amat and Pmat set in SNESSetJacobian() will be ignored
1224: -  mf_operator - use matrix-free only for the Amat used by SNESSetJacobian(), this means the user provided Pmat will continue to be used

1226:    Options Database:
1227: + -snes_mf - use matrix free for both the mat and pmat operator
1228: - -snes_mf_operator - use matrix free only for the mat operator

1230:    Level: intermediate

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

1234: .seealso:   SNESGetUseMatrixFree(), MatCreateSNESMF()
1235: @*/
1236: PetscErrorCode  SNESSetUseMatrixFree(SNES snes,PetscBool mf_operator,PetscBool mf)
1237: {
1242:   if (mf && !mf_operator) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"If using mf must also use mf_operator");
1243:   snes->mf          = mf;
1244:   snes->mf_operator = mf_operator;
1245:   return(0);
1246: }

1248: /*@
1249:    SNESGetUseMatrixFree - indicates if the SNES uses matrix free finite difference matrix vector products to apply
1250:                           the Jacobian.

1252:    Collective on SNES

1254:    Input Parameter:
1255: .  snes - SNES context

1257:    Output Parameters:
1258: +  mf - use matrix-free for both the Amat and Pmat used by SNESSetJacobian(), both the Amat and Pmat set in SNESSetJacobian() will be ignored
1259: -  mf_operator - use matrix-free only for the Amat used by SNESSetJacobian(), this means the user provided Pmat will continue to be used

1261:    Options Database:
1262: + -snes_mf - use matrix free for both the mat and pmat operator
1263: - -snes_mf_operator - use matrix free only for the mat operator

1265:    Level: intermediate

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

1269: .seealso:   SNESSetUseMatrixFree(), MatCreateSNESMF()
1270: @*/
1271: PetscErrorCode  SNESGetUseMatrixFree(SNES snes,PetscBool *mf_operator,PetscBool *mf)
1272: {
1275:   if (mf)          *mf          = snes->mf;
1276:   if (mf_operator) *mf_operator = snes->mf_operator;
1277:   return(0);
1278: }

1280: /*@
1281:    SNESGetIterationNumber - Gets the number of nonlinear iterations completed
1282:    at this time.

1284:    Not Collective

1286:    Input Parameter:
1287: .  snes - SNES context

1289:    Output Parameter:
1290: .  iter - iteration number

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

1295:    This is useful for using lagged Jacobians (where one does not recompute the
1296:    Jacobian at each SNES iteration). For example, the code
1297: .vb
1298:       SNESGetIterationNumber(snes,&it);
1299:       if (!(it % 2)) {
1300:         [compute Jacobian here]
1301:       }
1302: .ve
1303:    can be used in your ComputeJacobian() function to cause the Jacobian to be
1304:    recomputed every second SNES iteration.

1306:    After the SNES solve is complete this will return the number of nonlinear iterations used.

1308:    Level: intermediate

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

1312: .seealso:   SNESGetLinearSolveIterations()
1313: @*/
1314: PetscErrorCode  SNESGetIterationNumber(SNES snes,PetscInt *iter)
1315: {
1319:   *iter = snes->iter;
1320:   return(0);
1321: }

1323: /*@
1324:    SNESSetIterationNumber - Sets the current iteration number.

1326:    Not Collective

1328:    Input Parameter:
1329: .  snes - SNES context
1330: .  iter - iteration number

1332:    Level: developer

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

1336: .seealso:   SNESGetLinearSolveIterations()
1337: @*/
1338: PetscErrorCode  SNESSetIterationNumber(SNES snes,PetscInt iter)
1339: {

1344:   PetscObjectSAWsTakeAccess((PetscObject)snes);
1345:   snes->iter = iter;
1346:   PetscObjectSAWsGrantAccess((PetscObject)snes);
1347:   return(0);
1348: }

1350: /*@
1351:    SNESGetNonlinearStepFailures - Gets the number of unsuccessful steps
1352:    attempted by the nonlinear solver.

1354:    Not Collective

1356:    Input Parameter:
1357: .  snes - SNES context

1359:    Output Parameter:
1360: .  nfails - number of unsuccessful steps attempted

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

1365:    Level: intermediate

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

1369: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures(),
1370:           SNESSetMaxNonlinearStepFailures(), SNESGetMaxNonlinearStepFailures()
1371: @*/
1372: PetscErrorCode  SNESGetNonlinearStepFailures(SNES snes,PetscInt *nfails)
1373: {
1377:   *nfails = snes->numFailures;
1378:   return(0);
1379: }

1381: /*@
1382:    SNESSetMaxNonlinearStepFailures - Sets the maximum number of unsuccessful steps
1383:    attempted by the nonlinear solver before it gives up.

1385:    Not Collective

1387:    Input Parameters:
1388: +  snes     - SNES context
1389: -  maxFails - maximum of unsuccessful steps

1391:    Level: intermediate

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

1395: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures(),
1396:           SNESGetMaxNonlinearStepFailures(), SNESGetNonlinearStepFailures()
1397: @*/
1398: PetscErrorCode  SNESSetMaxNonlinearStepFailures(SNES snes, PetscInt maxFails)
1399: {
1402:   snes->maxFailures = maxFails;
1403:   return(0);
1404: }

1406: /*@
1407:    SNESGetMaxNonlinearStepFailures - Gets the maximum number of unsuccessful steps
1408:    attempted by the nonlinear solver before it gives up.

1410:    Not Collective

1412:    Input Parameter:
1413: .  snes     - SNES context

1415:    Output Parameter:
1416: .  maxFails - maximum of unsuccessful steps

1418:    Level: intermediate

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

1422: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures(),
1423:           SNESSetMaxNonlinearStepFailures(), SNESGetNonlinearStepFailures()

1425: @*/
1426: PetscErrorCode  SNESGetMaxNonlinearStepFailures(SNES snes, PetscInt *maxFails)
1427: {
1431:   *maxFails = snes->maxFailures;
1432:   return(0);
1433: }

1435: /*@
1436:    SNESGetNumberFunctionEvals - Gets the number of user provided function evaluations
1437:      done by SNES.

1439:    Not Collective

1441:    Input Parameter:
1442: .  snes     - SNES context

1444:    Output Parameter:
1445: .  nfuncs - number of evaluations

1447:    Level: intermediate

1449:    Notes:
1450:     Reset every time SNESSolve is called unless SNESSetCountersReset() is used.

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

1454: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures(), SNESSetCountersReset()
1455: @*/
1456: PetscErrorCode  SNESGetNumberFunctionEvals(SNES snes, PetscInt *nfuncs)
1457: {
1461:   *nfuncs = snes->nfuncs;
1462:   return(0);
1463: }

1465: /*@
1466:    SNESGetLinearSolveFailures - Gets the number of failed (non-converged)
1467:    linear solvers.

1469:    Not Collective

1471:    Input Parameter:
1472: .  snes - SNES context

1474:    Output Parameter:
1475: .  nfails - number of failed solves

1477:    Level: intermediate

1479:    Options Database Keys:
1480: . -snes_max_linear_solve_fail <num> - The number of failures before the solve is terminated

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

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

1487: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures()
1488: @*/
1489: PetscErrorCode  SNESGetLinearSolveFailures(SNES snes,PetscInt *nfails)
1490: {
1494:   *nfails = snes->numLinearSolveFailures;
1495:   return(0);
1496: }

1498: /*@
1499:    SNESSetMaxLinearSolveFailures - the number of failed linear solve attempts
1500:    allowed before SNES returns with a diverged reason of SNES_DIVERGED_LINEAR_SOLVE

1502:    Logically Collective on SNES

1504:    Input Parameters:
1505: +  snes     - SNES context
1506: -  maxFails - maximum allowed linear solve failures

1508:    Level: intermediate

1510:    Options Database Keys:
1511: . -snes_max_linear_solve_fail <num> - The number of failures before the solve is terminated

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

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

1518: .seealso: SNESGetLinearSolveFailures(), SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations()
1519: @*/
1520: PetscErrorCode  SNESSetMaxLinearSolveFailures(SNES snes, PetscInt maxFails)
1521: {
1525:   snes->maxLinearSolveFailures = maxFails;
1526:   return(0);
1527: }

1529: /*@
1530:    SNESGetMaxLinearSolveFailures - gets the maximum number of linear solve failures that
1531:      are allowed before SNES terminates

1533:    Not Collective

1535:    Input Parameter:
1536: .  snes     - SNES context

1538:    Output Parameter:
1539: .  maxFails - maximum of unsuccessful solves allowed

1541:    Level: intermediate

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

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

1548: .seealso: SNESGetLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(),
1549: @*/
1550: PetscErrorCode  SNESGetMaxLinearSolveFailures(SNES snes, PetscInt *maxFails)
1551: {
1555:   *maxFails = snes->maxLinearSolveFailures;
1556:   return(0);
1557: }

1559: /*@
1560:    SNESGetLinearSolveIterations - Gets the total number of linear iterations
1561:    used by the nonlinear solver.

1563:    Not Collective

1565:    Input Parameter:
1566: .  snes - SNES context

1568:    Output Parameter:
1569: .  lits - number of linear iterations

1571:    Notes:
1572:    This counter is reset to zero for each successive call to SNESSolve() unless SNESSetCountersReset() is used.

1574:    If the linear solver fails inside the SNESSolve() the iterations for that call to the linear solver are not included. If you wish to count them
1575:    then call KSPGetIterationNumber() after the failed solve.

1577:    Level: intermediate

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

1581: .seealso:  SNESGetIterationNumber(), SNESGetLinearSolveFailures(), SNESGetMaxLinearSolveFailures(), SNESSetCountersReset()
1582: @*/
1583: PetscErrorCode  SNESGetLinearSolveIterations(SNES snes,PetscInt *lits)
1584: {
1588:   *lits = snes->linear_its;
1589:   return(0);
1590: }

1592: /*@
1593:    SNESSetCountersReset - Sets whether or not the counters for linear iterations and function evaluations
1594:    are reset every time SNESSolve() is called.

1596:    Logically Collective on SNES

1598:    Input Parameter:
1599: +  snes - SNES context
1600: -  reset - whether to reset the counters or not

1602:    Notes:
1603:    This defaults to PETSC_TRUE

1605:    Level: developer

1607: .keywords: SNES, nonlinear, set, reset, number, linear, iterations

1609: .seealso:  SNESGetNumberFunctionEvals(), SNESGetLinearSolveIterations(), SNESGetNPC()
1610: @*/
1611: PetscErrorCode  SNESSetCountersReset(SNES snes,PetscBool reset)
1612: {
1616:   snes->counters_reset = reset;
1617:   return(0);
1618: }


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

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

1626:    Input Parameters:
1627: +  snes - the SNES context
1628: -  ksp - the KSP context

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

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

1637:    Level: developer

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

1641: .seealso: KSPGetPC(), SNESCreate(), KSPCreate(), SNESSetKSP()
1642: @*/
1643: PetscErrorCode  SNESSetKSP(SNES snes,KSP ksp)
1644: {

1651:   PetscObjectReference((PetscObject)ksp);
1652:   if (snes->ksp) {PetscObjectDereference((PetscObject)snes->ksp);}
1653:   snes->ksp = ksp;
1654:   return(0);
1655: }

1657: /* -----------------------------------------------------------*/
1658: /*@
1659:    SNESCreate - Creates a nonlinear solver context.

1661:    Collective on MPI_Comm

1663:    Input Parameters:
1664: .  comm - MPI communicator

1666:    Output Parameter:
1667: .  outsnes - the new SNES context

1669:    Options Database Keys:
1670: +   -snes_mf - Activates default matrix-free Jacobian-vector products,
1671:                and no preconditioning matrix
1672: .   -snes_mf_operator - Activates default matrix-free Jacobian-vector
1673:                products, and a user-provided preconditioning matrix
1674:                as set by SNESSetJacobian()
1675: -   -snes_fd - Uses (slow!) finite differences to compute Jacobian

1677:    Level: beginner

1679:    Developer Notes:
1680:     SNES always creates a KSP object even though many SNES methods do not use it. This is
1681:                     unfortunate and should be fixed at some point. The flag snes->usesksp indicates if the
1682:                     particular method does use KSP and regulates if the information about the KSP is printed
1683:                     in SNESView(). TSSetFromOptions() does call SNESSetFromOptions() which can lead to users being confused
1684:                     by help messages about meaningless SNES options.

1686:                     SNES always creates the snes->kspconvctx even though it is used by only one type. This should
1687:                     be fixed.

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

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

1693: @*/
1694: PetscErrorCode  SNESCreate(MPI_Comm comm,SNES *outsnes)
1695: {
1697:   SNES           snes;
1698:   SNESKSPEW      *kctx;

1702:   *outsnes = NULL;
1703:   SNESInitializePackage();

1705:   PetscHeaderCreate(snes,SNES_CLASSID,"SNES","Nonlinear solver","SNES",comm,SNESDestroy,SNESView);

1707:   snes->ops->converged    = SNESConvergedDefault;
1708:   snes->usesksp           = PETSC_TRUE;
1709:   snes->tolerancesset     = PETSC_FALSE;
1710:   snes->max_its           = 50;
1711:   snes->max_funcs         = 10000;
1712:   snes->norm              = 0.0;
1713:   snes->xnorm             = 0.0;
1714:   snes->ynorm             = 0.0;
1715:   snes->normschedule      = SNES_NORM_ALWAYS;
1716:   snes->functype          = SNES_FUNCTION_DEFAULT;
1717: #if defined(PETSC_USE_REAL_SINGLE)
1718:   snes->rtol              = 1.e-5;
1719: #else
1720:   snes->rtol              = 1.e-8;
1721: #endif
1722:   snes->ttol              = 0.0;
1723: #if defined(PETSC_USE_REAL_SINGLE)
1724:   snes->abstol            = 1.e-25;
1725: #else
1726:   snes->abstol            = 1.e-50;
1727: #endif
1728: #if defined(PETSC_USE_REAL_SINGLE)
1729:   snes->stol              = 1.e-5;
1730: #else
1731:   snes->stol              = 1.e-8;
1732: #endif
1733: #if defined(PETSC_USE_REAL_SINGLE)
1734:   snes->deltatol          = 1.e-6;
1735: #else
1736:   snes->deltatol          = 1.e-12;
1737: #endif
1738:   snes->divtol            = 1.e4;
1739:   snes->rnorm0            = 0;
1740:   snes->nfuncs            = 0;
1741:   snes->numFailures       = 0;
1742:   snes->maxFailures       = 1;
1743:   snes->linear_its        = 0;
1744:   snes->lagjacobian       = 1;
1745:   snes->jac_iter          = 0;
1746:   snes->lagjac_persist    = PETSC_FALSE;
1747:   snes->lagpreconditioner = 1;
1748:   snes->pre_iter          = 0;
1749:   snes->lagpre_persist    = PETSC_FALSE;
1750:   snes->numbermonitors    = 0;
1751:   snes->data              = 0;
1752:   snes->setupcalled       = PETSC_FALSE;
1753:   snes->ksp_ewconv        = PETSC_FALSE;
1754:   snes->nwork             = 0;
1755:   snes->work              = 0;
1756:   snes->nvwork            = 0;
1757:   snes->vwork             = 0;
1758:   snes->conv_hist_len     = 0;
1759:   snes->conv_hist_max     = 0;
1760:   snes->conv_hist         = NULL;
1761:   snes->conv_hist_its     = NULL;
1762:   snes->conv_hist_reset   = PETSC_TRUE;
1763:   snes->counters_reset    = PETSC_TRUE;
1764:   snes->vec_func_init_set = PETSC_FALSE;
1765:   snes->reason            = SNES_CONVERGED_ITERATING;
1766:   snes->npcside           = PC_RIGHT;
1767:   snes->setfromoptionscalled = 0;

1769:   snes->mf          = PETSC_FALSE;
1770:   snes->mf_operator = PETSC_FALSE;
1771:   snes->mf_version  = 1;

1773:   snes->numLinearSolveFailures = 0;
1774:   snes->maxLinearSolveFailures = 1;

1776:   snes->vizerotolerance = 1.e-8;
1777: #if defined(PETSC_USE_DEBUG)
1778:   snes->checkjacdomainerror = PETSC_TRUE;
1779: #else
1780:   snes->checkjacdomainerror = PETSC_FALSE;
1781: #endif

1783:   /* Set this to true if the implementation of SNESSolve_XXX does compute the residual at the final solution. */
1784:   snes->alwayscomputesfinalresidual = PETSC_FALSE;

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

1789:   snes->kspconvctx  = (void*)kctx;
1790:   kctx->version     = 2;
1791:   kctx->rtol_0      = .3; /* Eisenstat and Walker suggest rtol_0=.5, but
1792:                              this was too large for some test cases */
1793:   kctx->rtol_last   = 0.0;
1794:   kctx->rtol_max    = .9;
1795:   kctx->gamma       = 1.0;
1796:   kctx->alpha       = .5*(1.0 + PetscSqrtReal(5.0));
1797:   kctx->alpha2      = kctx->alpha;
1798:   kctx->threshold   = .1;
1799:   kctx->lresid_last = 0.0;
1800:   kctx->norm_last   = 0.0;

1802:   *outsnes = snes;
1803:   return(0);
1804: }

1806: /*MC
1807:     SNESFunction - Functional form used to convey the nonlinear function to be solved by SNES

1809:      Synopsis:
1810:      #include "petscsnes.h"
1811:      PetscErrorCode SNESFunction(SNES snes,Vec x,Vec f,void *ctx);

1813:      Input Parameters:
1814: +     snes - the SNES context
1815: .     x    - state at which to evaluate residual
1816: -     ctx     - optional user-defined function context, passed in with SNESSetFunction()

1818:      Output Parameter:
1819: .     f  - vector to put residual (function value)

1821:    Level: intermediate

1823: .seealso:   SNESSetFunction(), SNESGetFunction()
1824: M*/

1826: /*@C
1827:    SNESSetFunction - Sets the function evaluation routine and function
1828:    vector for use by the SNES routines in solving systems of nonlinear
1829:    equations.

1831:    Logically Collective on SNES

1833:    Input Parameters:
1834: +  snes - the SNES context
1835: .  r - vector to store function value
1836: .  f - function evaluation routine; see SNESFunction for calling sequence details
1837: -  ctx - [optional] user-defined context for private data for the
1838:          function evaluation routine (may be NULL)

1840:    Notes:
1841:    The Newton-like methods typically solve linear systems of the form
1842: $      f'(x) x = -f(x),
1843:    where f'(x) denotes the Jacobian matrix and f(x) is the function.

1845:    Level: beginner

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

1849: .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetPicard(), SNESFunction
1850: @*/
1851: PetscErrorCode  SNESSetFunction(SNES snes,Vec r,PetscErrorCode (*f)(SNES,Vec,Vec,void*),void *ctx)
1852: {
1854:   DM             dm;

1858:   if (r) {
1861:     PetscObjectReference((PetscObject)r);
1862:     VecDestroy(&snes->vec_func);

1864:     snes->vec_func = r;
1865:   }
1866:   SNESGetDM(snes,&dm);
1867:   DMSNESSetFunction(dm,f,ctx);
1868:   return(0);
1869: }


1872: /*@C
1873:    SNESSetInitialFunction - Sets the function vector to be used as the
1874:    function norm at the initialization of the method.  In some
1875:    instances, the user has precomputed the function before calling
1876:    SNESSolve.  This function allows one to avoid a redundant call
1877:    to SNESComputeFunction in that case.

1879:    Logically Collective on SNES

1881:    Input Parameters:
1882: +  snes - the SNES context
1883: -  f - vector to store function value

1885:    Notes:
1886:    This should not be modified during the solution procedure.

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

1890:    Level: developer

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

1894: .seealso: SNESSetFunction(), SNESComputeFunction(), SNESSetInitialFunctionNorm()
1895: @*/
1896: PetscErrorCode  SNESSetInitialFunction(SNES snes, Vec f)
1897: {
1899:   Vec            vec_func;

1905:   if (snes->npcside== PC_LEFT && snes->functype == SNES_FUNCTION_PRECONDITIONED) {
1906:     snes->vec_func_init_set = PETSC_FALSE;
1907:     return(0);
1908:   }
1909:   SNESGetFunction(snes,&vec_func,NULL,NULL);
1910:   VecCopy(f, vec_func);

1912:   snes->vec_func_init_set = PETSC_TRUE;
1913:   return(0);
1914: }

1916: /*@
1917:    SNESSetNormSchedule - Sets the SNESNormSchedule used in covergence and monitoring
1918:    of the SNES method.

1920:    Logically Collective on SNES

1922:    Input Parameters:
1923: +  snes - the SNES context
1924: -  normschedule - the frequency of norm computation

1926:    Options Database Key:
1927: .  -snes_norm_schedule <none, always, initialonly, finalonly, initalfinalonly>

1929:    Notes:
1930:    Only certain SNES methods support certain SNESNormSchedules.  Most require evaluation
1931:    of the nonlinear function and the taking of its norm at every iteration to
1932:    even ensure convergence at all.  However, methods such as custom Gauss-Seidel methods
1933:    (SNESNGS) and the like do not require the norm of the function to be computed, and therfore
1934:    may either be monitored for convergence or not.  As these are often used as nonlinear
1935:    preconditioners, monitoring the norm of their error is not a useful enterprise within
1936:    their solution.

1938:    Level: developer

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

1942: .seealso: SNESGetNormSchedule(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormSchedule
1943: @*/
1944: PetscErrorCode  SNESSetNormSchedule(SNES snes, SNESNormSchedule normschedule)
1945: {
1948:   snes->normschedule = normschedule;
1949:   return(0);
1950: }


1953: /*@
1954:    SNESGetNormSchedule - Gets the SNESNormSchedule used in covergence and monitoring
1955:    of the SNES method.

1957:    Logically Collective on SNES

1959:    Input Parameters:
1960: +  snes - the SNES context
1961: -  normschedule - the type of the norm used

1963:    Level: advanced

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

1967: .seealso: SNESSetNormSchedule(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormSchedule
1968: @*/
1969: PetscErrorCode  SNESGetNormSchedule(SNES snes, SNESNormSchedule *normschedule)
1970: {
1973:   *normschedule = snes->normschedule;
1974:   return(0);
1975: }


1978: /*@
1979:   SNESSetFunctionNorm - Sets the last computed residual norm.

1981:   Logically Collective on SNES

1983:   Input Parameters:
1984: + snes - the SNES context

1986: - normschedule - the frequency of norm computation

1988:   Level: developer

1990: .keywords: SNES, nonlinear, set, function, norm, type
1991: .seealso: SNESGetNormSchedule(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormSchedule
1992: @*/
1993: PetscErrorCode SNESSetFunctionNorm(SNES snes, PetscReal norm)
1994: {
1997:   snes->norm = norm;
1998:   return(0);
1999: }

2001: /*@
2002:   SNESGetFunctionNorm - Gets the last computed norm of the residual

2004:   Not Collective

2006:   Input Parameter:
2007: . snes - the SNES context

2009:   Output Parameter:
2010: . norm - the last computed residual norm

2012:   Level: developer

2014: .keywords: SNES, nonlinear, set, function, norm, type
2015: .seealso: SNESSetNormSchedule(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormSchedule
2016: @*/
2017: PetscErrorCode SNESGetFunctionNorm(SNES snes, PetscReal *norm)
2018: {
2022:   *norm = snes->norm;
2023:   return(0);
2024: }

2026: /*@
2027:   SNESGetUpdateNorm - Gets the last computed norm of the Newton update

2029:   Not Collective

2031:   Input Parameter:
2032: . snes - the SNES context

2034:   Output Parameter:
2035: . ynorm - the last computed update norm

2037:   Level: developer

2039: .keywords: SNES, nonlinear, set, function, norm, type
2040: .seealso: SNESSetNormSchedule(), SNESComputeFunction(), SNESGetFunctionNorm()
2041: @*/
2042: PetscErrorCode SNESGetUpdateNorm(SNES snes, PetscReal *ynorm)
2043: {
2047:   *ynorm = snes->ynorm;
2048:   return(0);
2049: }

2051: /*@
2052:   SNESGetSolutionNorm - Gets the last computed norm of the solution

2054:   Not Collective

2056:   Input Parameter:
2057: . snes - the SNES context

2059:   Output Parameter:
2060: . xnorm - the last computed solution norm

2062:   Level: developer

2064: .keywords: SNES, nonlinear, set, function, norm, type
2065: .seealso: SNESSetNormSchedule(), SNESComputeFunction(), SNESGetFunctionNorm(), SNESGetUpdateNorm()
2066: @*/
2067: PetscErrorCode SNESGetSolutionNorm(SNES snes, PetscReal *xnorm)
2068: {
2072:   *xnorm = snes->xnorm;
2073:   return(0);
2074: }

2076: /*@C
2077:    SNESSetFunctionType - Sets the SNESNormSchedule used in covergence and monitoring
2078:    of the SNES method.

2080:    Logically Collective on SNES

2082:    Input Parameters:
2083: +  snes - the SNES context
2084: -  normschedule - the frequency of norm computation

2086:    Notes:
2087:    Only certain SNES methods support certain SNESNormSchedules.  Most require evaluation
2088:    of the nonlinear function and the taking of its norm at every iteration to
2089:    even ensure convergence at all.  However, methods such as custom Gauss-Seidel methods
2090:    (SNESNGS) and the like do not require the norm of the function to be computed, and therfore
2091:    may either be monitored for convergence or not.  As these are often used as nonlinear
2092:    preconditioners, monitoring the norm of their error is not a useful enterprise within
2093:    their solution.

2095:    Level: developer

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

2099: .seealso: SNESGetNormSchedule(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormSchedule
2100: @*/
2101: PetscErrorCode  SNESSetFunctionType(SNES snes, SNESFunctionType type)
2102: {
2105:   snes->functype = type;
2106:   return(0);
2107: }


2110: /*@C
2111:    SNESGetFunctionType - Gets the SNESNormSchedule used in covergence and monitoring
2112:    of the SNES method.

2114:    Logically Collective on SNES

2116:    Input Parameters:
2117: +  snes - the SNES context
2118: -  normschedule - the type of the norm used

2120:    Level: advanced

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

2124: .seealso: SNESSetNormSchedule(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormSchedule
2125: @*/
2126: PetscErrorCode  SNESGetFunctionType(SNES snes, SNESFunctionType *type)
2127: {
2130:   *type = snes->functype;
2131:   return(0);
2132: }

2134: /*MC
2135:     SNESNGSFunction - function used to convey a Gauss-Seidel sweep on the nonlinear function

2137:      Synopsis:
2138:      #include <petscsnes.h>
2139: $    SNESNGSFunction(SNES snes,Vec x,Vec b,void *ctx);

2141: +  X   - solution vector
2142: .  B   - RHS vector
2143: -  ctx - optional user-defined Gauss-Seidel context

2145:    Level: intermediate

2147: .seealso:   SNESSetNGS(), SNESGetNGS()
2148: M*/

2150: /*@C
2151:    SNESSetNGS - Sets the user nonlinear Gauss-Seidel routine for
2152:    use with composed nonlinear solvers.

2154:    Input Parameters:
2155: +  snes   - the SNES context
2156: .  f - function evaluation routine to apply Gauss-Seidel see SNESNGSFunction
2157: -  ctx    - [optional] user-defined context for private data for the
2158:             smoother evaluation routine (may be NULL)

2160:    Notes:
2161:    The NGS routines are used by the composed nonlinear solver to generate
2162:     a problem appropriate update to the solution, particularly FAS.

2164:    Level: intermediate

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

2168: .seealso: SNESGetFunction(), SNESComputeNGS()
2169: @*/
2170: PetscErrorCode SNESSetNGS(SNES snes,PetscErrorCode (*f)(SNES,Vec,Vec,void*),void *ctx)
2171: {
2173:   DM             dm;

2177:   SNESGetDM(snes,&dm);
2178:   DMSNESSetNGS(dm,f,ctx);
2179:   return(0);
2180: }

2182: PetscErrorCode SNESPicardComputeFunction(SNES snes,Vec x,Vec f,void *ctx)
2183: {
2185:   DM             dm;
2186:   DMSNES         sdm;

2189:   SNESGetDM(snes,&dm);
2190:   DMGetDMSNES(dm,&sdm);
2191:   if (!sdm->ops->computepfunction) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetPicard() to provide Picard function.");
2192:   if (!sdm->ops->computepjacobian) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetPicard() to provide Picard Jacobian.");
2193:   /*  A(x)*x - b(x) */
2194:   PetscStackPush("SNES Picard user function");
2195:   (*sdm->ops->computepfunction)(snes,x,f,sdm->pctx);
2196:   PetscStackPop;
2197:   PetscStackPush("SNES Picard user Jacobian");
2198:   (*sdm->ops->computepjacobian)(snes,x,snes->jacobian,snes->jacobian_pre,sdm->pctx);
2199:   PetscStackPop;
2200:   VecScale(f,-1.0);
2201:   MatMultAdd(snes->jacobian,x,f,f);
2202:   return(0);
2203: }

2205: PetscErrorCode SNESPicardComputeJacobian(SNES snes,Vec x1,Mat J,Mat B,void *ctx)
2206: {
2208:   /* the jacobian matrix should be pre-filled in SNESPicardComputeFunction */
2209:   return(0);
2210: }

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

2215:    Logically Collective on SNES

2217:    Input Parameters:
2218: +  snes - the SNES context
2219: .  r - vector to store function value
2220: .  b - function evaluation routine
2221: .  Amat - matrix with which A(x) x - b(x) is to be computed
2222: .  Pmat - matrix from which preconditioner is computed (usually the same as Amat)
2223: .  J  - function to compute matrix value, see SNESJacobianFunction for details on its calling sequence
2224: -  ctx - [optional] user-defined context for private data for the
2225:          function evaluation routine (may be NULL)

2227:    Notes:
2228:     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
2229:     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.

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

2233: $     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}
2234: $     Note that when an exact solver is used this corresponds to the "classic" Picard A(x^{n}) x^{n+1} = b(x^{n}) iteration.

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

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

2241:    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
2242:    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
2243:    different please contact us at petsc-dev@mcs.anl.gov and we'll have an entirely new argument :-).

2245:    Level: intermediate

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

2249: .seealso: SNESGetFunction(), SNESSetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESGetPicard(), SNESLineSearchPreCheckPicard(), SNESJacobianFunction
2250: @*/
2251: PetscErrorCode  SNESSetPicard(SNES snes,Vec r,PetscErrorCode (*b)(SNES,Vec,Vec,void*),Mat Amat, Mat Pmat, PetscErrorCode (*J)(SNES,Vec,Mat,Mat,void*),void *ctx)
2252: {
2254:   DM             dm;

2258:   SNESGetDM(snes, &dm);
2259:   DMSNESSetPicard(dm,b,J,ctx);
2260:   SNESSetFunction(snes,r,SNESPicardComputeFunction,ctx);
2261:   SNESSetJacobian(snes,Amat,Pmat,SNESPicardComputeJacobian,ctx);
2262:   return(0);
2263: }

2265: /*@C
2266:    SNESGetPicard - Returns the context for the Picard iteration

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

2270:    Input Parameter:
2271: .  snes - the SNES context

2273:    Output Parameter:
2274: +  r - the function (or NULL)
2275: .  f - the function (or NULL); see SNESFunction for calling sequence details
2276: .  Amat - the matrix used to defined the operation A(x) x - b(x) (or NULL)
2277: .  Pmat  - the matrix from which the preconditioner will be constructed (or NULL)
2278: .  J - the function for matrix evaluation (or NULL); see SNESJacobianFunction for calling sequence details
2279: -  ctx - the function context (or NULL)

2281:    Level: advanced

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

2285: .seealso: SNESSetPicard(), SNESGetFunction(), SNESGetJacobian(), SNESGetDM(), SNESFunction, SNESJacobianFunction
2286: @*/
2287: PetscErrorCode  SNESGetPicard(SNES snes,Vec *r,PetscErrorCode (**f)(SNES,Vec,Vec,void*),Mat *Amat, Mat *Pmat, PetscErrorCode (**J)(SNES,Vec,Mat,Mat,void*),void **ctx)
2288: {
2290:   DM             dm;

2294:   SNESGetFunction(snes,r,NULL,NULL);
2295:   SNESGetJacobian(snes,Amat,Pmat,NULL,NULL);
2296:   SNESGetDM(snes,&dm);
2297:   DMSNESGetPicard(dm,f,J,ctx);
2298:   return(0);
2299: }

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

2304:    Logically Collective on SNES

2306:    Input Parameters:
2307: +  snes - the SNES context
2308: .  func - function evaluation routine
2309: -  ctx - [optional] user-defined context for private data for the
2310:          function evaluation routine (may be NULL)

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

2315: .  f - function vector
2316: -  ctx - optional user-defined function context

2318:    Level: intermediate

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

2322: .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian()
2323: @*/
2324: PetscErrorCode  SNESSetComputeInitialGuess(SNES snes,PetscErrorCode (*func)(SNES,Vec,void*),void *ctx)
2325: {
2328:   if (func) snes->ops->computeinitialguess = func;
2329:   if (ctx)  snes->initialguessP            = ctx;
2330:   return(0);
2331: }

2333: /* --------------------------------------------------------------- */
2334: /*@C
2335:    SNESGetRhs - Gets the vector for solving F(x) = rhs. If rhs is not set
2336:    it assumes a zero right hand side.

2338:    Logically Collective on SNES

2340:    Input Parameter:
2341: .  snes - the SNES context

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

2346:    Level: intermediate

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

2350: .seealso: SNESGetSolution(), SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction()
2351: @*/
2352: PetscErrorCode  SNESGetRhs(SNES snes,Vec *rhs)
2353: {
2357:   *rhs = snes->vec_rhs;
2358:   return(0);
2359: }

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

2364:    Collective on SNES

2366:    Input Parameters:
2367: +  snes - the SNES context
2368: -  x - input vector

2370:    Output Parameter:
2371: .  y - function vector, as set by SNESSetFunction()

2373:    Notes:
2374:    SNESComputeFunction() is typically used within nonlinear solvers
2375:    implementations, so most users would not generally call this routine
2376:    themselves.

2378:    Level: developer

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

2382: .seealso: SNESSetFunction(), SNESGetFunction()
2383: @*/
2384: PetscErrorCode  SNESComputeFunction(SNES snes,Vec x,Vec y)
2385: {
2387:   DM             dm;
2388:   DMSNES         sdm;

2396:   VecValidValues(x,2,PETSC_TRUE);

2398:   SNESGetDM(snes,&dm);
2399:   DMGetDMSNES(dm,&sdm);
2400:   if (sdm->ops->computefunction) {
2401:     if (sdm->ops->computefunction != SNESObjectiveComputeFunctionDefaultFD) {
2402:       PetscLogEventBegin(SNES_FunctionEval,snes,x,y,0);
2403:     }
2404:     VecLockReadPush(x);
2405:     PetscStackPush("SNES user function");
2406:     (*sdm->ops->computefunction)(snes,x,y,sdm->functionctx);
2407:     PetscStackPop;
2408:     VecLockReadPop(x);
2409:     if (sdm->ops->computefunction != SNESObjectiveComputeFunctionDefaultFD) {
2410:       PetscLogEventEnd(SNES_FunctionEval,snes,x,y,0);
2411:     }
2412:   } else if (snes->vec_rhs) {
2413:     MatMult(snes->jacobian, x, y);
2414:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetFunction() or SNESSetDM() before SNESComputeFunction(), likely called from SNESSolve().");
2415:   if (snes->vec_rhs) {
2416:     VecAXPY(y,-1.0,snes->vec_rhs);
2417:   }
2418:   snes->nfuncs++;
2419:   /*
2420:      domainerror might not be set on all processes; so we tag vector locally with Inf and the next inner product or norm will
2421:      propagate the value to all processes
2422:   */
2423:   if (snes->domainerror) {
2424:     VecSetInf(y);
2425:   }
2426:   return(0);
2427: }

2429: /*@
2430:    SNESComputeNGS - Calls the Gauss-Seidel function that has been set with  SNESSetNGS().

2432:    Collective on SNES

2434:    Input Parameters:
2435: +  snes - the SNES context
2436: .  x - input vector
2437: -  b - rhs vector

2439:    Output Parameter:
2440: .  x - new solution vector

2442:    Notes:
2443:    SNESComputeNGS() is typically used within composed nonlinear solver
2444:    implementations, so most users would not generally call this routine
2445:    themselves.

2447:    Level: developer

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

2451: .seealso: SNESSetNGS(), SNESComputeFunction()
2452: @*/
2453: PetscErrorCode  SNESComputeNGS(SNES snes,Vec b,Vec x)
2454: {
2456:   DM             dm;
2457:   DMSNES         sdm;

2465:   if (b) {VecValidValues(b,2,PETSC_TRUE);}
2466:   PetscLogEventBegin(SNES_NGSEval,snes,x,b,0);
2467:   SNESGetDM(snes,&dm);
2468:   DMGetDMSNES(dm,&sdm);
2469:   if (sdm->ops->computegs) {
2470:     if (b) {VecLockReadPush(b);}
2471:     PetscStackPush("SNES user NGS");
2472:     (*sdm->ops->computegs)(snes,x,b,sdm->gsctx);
2473:     PetscStackPop;
2474:     if (b) {VecLockReadPop(b);}
2475:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetNGS() before SNESComputeNGS(), likely called from SNESSolve().");
2476:   PetscLogEventEnd(SNES_NGSEval,snes,x,b,0);
2477:   return(0);
2478: }

2480: PetscErrorCode SNESTestJacobian(SNES snes)
2481: {
2482:   Mat               A,B,C,D,jacobian;
2483:   Vec               x = snes->vec_sol,f = snes->vec_func;
2484:   PetscErrorCode    ierr;
2485:   PetscReal         nrm,gnorm;
2486:   PetscReal         threshold = 1.e-5;
2487:   PetscInt          m,n,M,N;
2488:   void              *functx;
2489:   PetscBool         complete_print = PETSC_FALSE,threshold_print = PETSC_FALSE,test = PETSC_FALSE,flg;
2490:   PetscViewer       viewer,mviewer;
2491:   MPI_Comm          comm;
2492:   PetscInt          tabs;
2493:   static PetscBool  directionsprinted = PETSC_FALSE;
2494:   PetscViewerFormat format;

2497:   PetscObjectOptionsBegin((PetscObject)snes);
2498:   PetscOptionsName("-snes_test_jacobian","Compare hand-coded and finite difference Jacobians","None",&test);
2499:   PetscOptionsReal("-snes_test_jacobian", "Threshold for element difference between hand-coded and finite difference being meaningful", "None", threshold, &threshold,NULL);
2500:   PetscOptionsViewer("-snes_test_jacobian_view","View difference between hand-coded and finite difference Jacobians element entries","None",&mviewer,&format,&complete_print);
2501:   if (!complete_print) {
2502:     PetscOptionsViewer("-snes_test_jacobian_display","Display difference between hand-coded and finite difference Jacobians","None",&mviewer,&format,&complete_print);
2503:   }
2504:   /* for compatibility with PETSc 3.9 and older. */
2505:   PetscOptionsReal("-snes_test_jacobian_display_threshold", "Display difference between hand-coded and finite difference Jacobians which exceed input threshold", "None", threshold, &threshold, &threshold_print);
2506:   PetscOptionsEnd();
2507:   if (!test) return(0);

2509:   PetscObjectGetComm((PetscObject)snes,&comm);
2510:   PetscViewerASCIIGetStdout(comm,&viewer);
2511:   PetscViewerASCIIGetTab(viewer, &tabs);
2512:   PetscViewerASCIISetTab(viewer, ((PetscObject)snes)->tablevel);
2513:   PetscViewerASCIIPrintf(viewer,"  ---------- Testing Jacobian -------------\n");
2514:   if (!complete_print && !directionsprinted) {
2515:     PetscViewerASCIIPrintf(viewer,"  Run with -snes_test_jacobian_view and optionally -snes_test_jacobian <threshold> to show difference\n");
2516:     PetscViewerASCIIPrintf(viewer,"    of hand-coded and finite difference Jacobian entries greater than <threshold>.\n");
2517:   }
2518:   if (!directionsprinted) {
2519:     PetscViewerASCIIPrintf(viewer,"  Testing hand-coded Jacobian, if (for double precision runs) ||J - Jfd||_F/||J||_F is\n");
2520:     PetscViewerASCIIPrintf(viewer,"    O(1.e-8), the hand-coded Jacobian is probably correct.\n");
2521:     directionsprinted = PETSC_TRUE;
2522:   }
2523:   if (complete_print) {
2524:     PetscViewerPushFormat(mviewer,format);
2525:   }

2527:   /* evaluate the function at this point because SNESComputeJacobianDefault() assumes that the function has been evaluated and put into snes->vec_func */
2528:   SNESComputeFunction(snes,x,f);

2530:   PetscObjectTypeCompare((PetscObject)snes->jacobian,MATMFFD,&flg);
2531:   if (!flg) jacobian = snes->jacobian;
2532:   else jacobian = snes->jacobian_pre;

2534:   while (jacobian) {
2535:     PetscObjectTypeCompareAny((PetscObject)jacobian,&flg,MATSEQAIJ,MATMPIAIJ,MATSEQDENSE,MATMPIDENSE,MATSEQBAIJ,MATMPIBAIJ,MATSEQSBAIJ,MATMPIBAIJ,"");
2536:     if (flg) {
2537:       A    = jacobian;
2538:       PetscObjectReference((PetscObject)A);
2539:     } else {
2540:       MatComputeExplicitOperator(jacobian,&A);
2541:     }

2543:     MatCreate(PetscObjectComm((PetscObject)A),&B);
2544:     MatGetSize(A,&M,&N);
2545:     MatGetLocalSize(A,&m,&n);
2546:     MatSetSizes(B,m,n,M,N);
2547:     MatSetType(B,((PetscObject)A)->type_name);
2548:     MatSetUp(B);
2549:     MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);

2551:     SNESGetFunction(snes,NULL,NULL,&functx);
2552:     SNESComputeJacobianDefault(snes,x,B,B,functx);

2554:     MatDuplicate(B,MAT_COPY_VALUES,&D);
2555:     MatAYPX(D,-1.0,A,DIFFERENT_NONZERO_PATTERN);
2556:     MatNorm(D,NORM_FROBENIUS,&nrm);
2557:     MatNorm(A,NORM_FROBENIUS,&gnorm);
2558:     MatDestroy(&D);
2559:     if (!gnorm) gnorm = 1; /* just in case */
2560:     PetscViewerASCIIPrintf(viewer,"  ||J - Jfd||_F/||J||_F = %g, ||J - Jfd||_F = %g\n",(double)(nrm/gnorm),(double)nrm);

2562:     if (complete_print) {
2563:       PetscViewerASCIIPrintf(viewer,"  Hand-coded Jacobian ----------\n");
2564:       MatView(jacobian,mviewer);
2565:       PetscViewerASCIIPrintf(viewer,"  Finite difference Jacobian ----------\n");
2566:       MatView(B,mviewer);
2567:     }

2569:     if (threshold_print || complete_print) {
2570:       PetscInt          Istart, Iend, *ccols, bncols, cncols, j, row;
2571:       PetscScalar       *cvals;
2572:       const PetscInt    *bcols;
2573:       const PetscScalar *bvals;

2575:       MatAYPX(B,-1.0,A,DIFFERENT_NONZERO_PATTERN);
2576:       MatCreate(PetscObjectComm((PetscObject)A),&C);
2577:       MatSetSizes(C,m,n,M,N);
2578:       MatSetType(C,((PetscObject)A)->type_name);
2579:       MatSetUp(C);
2580:       MatSetOption(C,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
2581:       MatGetOwnershipRange(B,&Istart,&Iend);

2583:       for (row = Istart; row < Iend; row++) {
2584:         MatGetRow(B,row,&bncols,&bcols,&bvals);
2585:         PetscMalloc2(bncols,&ccols,bncols,&cvals);
2586:         for (j = 0, cncols = 0; j < bncols; j++) {
2587:           if (PetscAbsScalar(bvals[j]) > threshold) {
2588:             ccols[cncols] = bcols[j];
2589:             cvals[cncols] = bvals[j];
2590:             cncols += 1;
2591:           }
2592:         }
2593:         if (cncols) {
2594:           MatSetValues(C,1,&row,cncols,ccols,cvals,INSERT_VALUES);
2595:         }
2596:         MatRestoreRow(B,row,&bncols,&bcols,&bvals);
2597:         PetscFree2(ccols,cvals);
2598:       }
2599:       MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
2600:       MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
2601:       PetscViewerASCIIPrintf(viewer,"  Hand-coded minus finite-difference Jacobian with tolerance %g ----------\n",(double)threshold);
2602:       MatView(C,complete_print ? mviewer : viewer);
2603:       MatDestroy(&C);
2604:     }
2605:     MatDestroy(&A);
2606:     MatDestroy(&B);

2608:     if (jacobian != snes->jacobian_pre) {
2609:       jacobian = snes->jacobian_pre;
2610:       PetscViewerASCIIPrintf(viewer,"  ---------- Testing Jacobian for preconditioner -------------\n");
2611:     }
2612:     else jacobian = NULL;
2613:   }
2614:   if (complete_print) {
2615:     PetscViewerPopFormat(mviewer);
2616:   }
2617:   if (mviewer) { PetscViewerDestroy(&mviewer); }
2618:   PetscViewerASCIISetTab(viewer,tabs);
2619:   return(0);
2620: }

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

2625:    Collective on SNES and Mat

2627:    Input Parameters:
2628: +  snes - the SNES context
2629: -  x - input vector

2631:    Output Parameters:
2632: +  A - Jacobian matrix
2633: -  B - optional preconditioning matrix

2635:   Options Database Keys:
2636: +    -snes_lag_preconditioner <lag>
2637: .    -snes_lag_jacobian <lag>
2638: .    -snes_test_jacobian - compare the user provided Jacobian with one compute via finite differences to check for errors
2639: .    -snes_test_jacobian_display - display the user provided Jacobian, the finite difference Jacobian and the difference between them to help users detect the location of errors in the user provided Jacobian
2640: .    -snes_test_jacobian_display_threshold <numerical value>  - display entries in the difference between the user provided Jacobian and finite difference Jacobian that are greater than a certain value to help users detect errors
2641: .    -snes_compare_explicit - Compare the computed Jacobian to the finite difference Jacobian and output the differences
2642: .    -snes_compare_explicit_draw  - Compare the computed Jacobian to the finite difference Jacobian and draw the result
2643: .    -snes_compare_explicit_contour  - Compare the computed Jacobian to the finite difference Jacobian and draw a contour plot with the result
2644: .    -snes_compare_operator  - Make the comparison options above use the operator instead of the preconditioning matrix
2645: .    -snes_compare_coloring - Compute the finite difference Jacobian using coloring and display norms of difference
2646: .    -snes_compare_coloring_display - Compute the finite differece Jacobian using coloring and display verbose differences
2647: .    -snes_compare_coloring_threshold - Display only those matrix entries that differ by more than a given threshold
2648: .    -snes_compare_coloring_threshold_atol - Absolute tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold
2649: .    -snes_compare_coloring_threshold_rtol - Relative tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold
2650: .    -snes_compare_coloring_draw - Compute the finite differece Jacobian using coloring and draw differences
2651: -    -snes_compare_coloring_draw_contour - Compute the finite differece Jacobian using coloring and show contours of matrices and differences


2654:    Notes:
2655:    Most users should not need to explicitly call this routine, as it
2656:    is used internally within the nonlinear solvers.

2658:    Developer Notes:
2659:     This has duplicative ways of checking the accuracy of the user provided Jacobian (see the options above). This is for historical reasons, the routine SNESTestJacobian() use to used 
2660:       for with the SNESType of test that has been removed.

2662:    Level: developer

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

2666: .seealso:  SNESSetJacobian(), KSPSetOperators(), MatStructure, SNESSetLagPreconditioner(), SNESSetLagJacobian()
2667: @*/
2668: PetscErrorCode  SNESComputeJacobian(SNES snes,Vec X,Mat A,Mat B)
2669: {
2671:   PetscBool      flag;
2672:   DM             dm;
2673:   DMSNES         sdm;
2674:   KSP            ksp;

2680:   VecValidValues(X,2,PETSC_TRUE);
2681:   SNESGetDM(snes,&dm);
2682:   DMGetDMSNES(dm,&sdm);

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

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

2688:   if (snes->lagjacobian == -2) {
2689:     snes->lagjacobian = -1;

2691:     PetscInfo(snes,"Recomputing Jacobian/preconditioner because lag is -2 (means compute Jacobian, but then never again) \n");
2692:   } else if (snes->lagjacobian == -1) {
2693:     PetscInfo(snes,"Reusing Jacobian/preconditioner because lag is -1\n");
2694:     PetscObjectTypeCompare((PetscObject)A,MATMFFD,&flag);
2695:     if (flag) {
2696:       MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
2697:       MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
2698:     }
2699:     return(0);
2700:   } else if (snes->lagjacobian > 1 && (snes->iter + snes->jac_iter) % snes->lagjacobian) {
2701:     PetscInfo2(snes,"Reusing Jacobian/preconditioner because lag is %D and SNES iteration is %D\n",snes->lagjacobian,snes->iter);
2702:     PetscObjectTypeCompare((PetscObject)A,MATMFFD,&flag);
2703:     if (flag) {
2704:       MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
2705:       MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
2706:     }
2707:     return(0);
2708:   }
2709:   if (snes->npc && snes->npcside== PC_LEFT) {
2710:       MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
2711:       MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
2712:       return(0);
2713:   }

2715:   PetscLogEventBegin(SNES_JacobianEval,snes,X,A,B);
2716:   VecLockReadPush(X);
2717:   PetscStackPush("SNES user Jacobian function");
2718:   (*sdm->ops->computejacobian)(snes,X,A,B,sdm->jacobianctx);
2719:   PetscStackPop;
2720:   VecLockReadPop(X);
2721:   PetscLogEventEnd(SNES_JacobianEval,snes,X,A,B);

2723:   /* the next line ensures that snes->ksp exists */
2724:   SNESGetKSP(snes,&ksp);
2725:   if (snes->lagpreconditioner == -2) {
2726:     PetscInfo(snes,"Rebuilding preconditioner exactly once since lag is -2\n");
2727:     KSPSetReusePreconditioner(snes->ksp,PETSC_FALSE);
2728:     snes->lagpreconditioner = -1;
2729:   } else if (snes->lagpreconditioner == -1) {
2730:     PetscInfo(snes,"Reusing preconditioner because lag is -1\n");
2731:     KSPSetReusePreconditioner(snes->ksp,PETSC_TRUE);
2732:   } else if (snes->lagpreconditioner > 1 && (snes->iter + snes->pre_iter) % snes->lagpreconditioner) {
2733:     PetscInfo2(snes,"Reusing preconditioner because lag is %D and SNES iteration is %D\n",snes->lagpreconditioner,snes->iter);
2734:     KSPSetReusePreconditioner(snes->ksp,PETSC_TRUE);
2735:   } else {
2736:     PetscInfo(snes,"Rebuilding preconditioner\n");
2737:     KSPSetReusePreconditioner(snes->ksp,PETSC_FALSE);
2738:   }

2740:   SNESTestJacobian(snes);
2741:   /* make sure user returned a correct Jacobian and preconditioner */
2744:   {
2745:     PetscBool flag = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_operator = PETSC_FALSE;
2746:     PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject) snes)->options,((PetscObject)snes)->prefix,"-snes_compare_explicit",NULL,NULL,&flag);
2747:     PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject) snes)->options,((PetscObject)snes)->prefix,"-snes_compare_explicit_draw",NULL,NULL,&flag_draw);
2748:     PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject) snes)->options,((PetscObject)snes)->prefix,"-snes_compare_explicit_draw_contour",NULL,NULL,&flag_contour);
2749:     PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject) snes)->options,((PetscObject)snes)->prefix,"-snes_compare_operator",NULL,NULL,&flag_operator);
2750:     if (flag || flag_draw || flag_contour) {
2751:       Mat          Bexp_mine = NULL,Bexp,FDexp;
2752:       PetscViewer  vdraw,vstdout;
2753:       PetscBool    flg;
2754:       if (flag_operator) {
2755:         MatComputeExplicitOperator(A,&Bexp_mine);
2756:         Bexp = Bexp_mine;
2757:       } else {
2758:         /* See if the preconditioning matrix can be viewed and added directly */
2759:         PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQAIJ,MATMPIAIJ,MATSEQDENSE,MATMPIDENSE,MATSEQBAIJ,MATMPIBAIJ,MATSEQSBAIJ,MATMPIBAIJ,"");
2760:         if (flg) Bexp = B;
2761:         else {
2762:           /* If the "preconditioning" matrix is itself MATSHELL or some other type without direct support */
2763:           MatComputeExplicitOperator(B,&Bexp_mine);
2764:           Bexp = Bexp_mine;
2765:         }
2766:       }
2767:       MatConvert(Bexp,MATSAME,MAT_INITIAL_MATRIX,&FDexp);
2768:       SNESComputeJacobianDefault(snes,X,FDexp,FDexp,NULL);
2769:       PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&vstdout);
2770:       if (flag_draw || flag_contour) {
2771:         PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes),0,"Explicit Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);
2772:         if (flag_contour) {PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);}
2773:       } else vdraw = NULL;
2774:       PetscViewerASCIIPrintf(vstdout,"Explicit %s\n",flag_operator ? "Jacobian" : "preconditioning Jacobian");
2775:       if (flag) {MatView(Bexp,vstdout);}
2776:       if (vdraw) {MatView(Bexp,vdraw);}
2777:       PetscViewerASCIIPrintf(vstdout,"Finite difference Jacobian\n");
2778:       if (flag) {MatView(FDexp,vstdout);}
2779:       if (vdraw) {MatView(FDexp,vdraw);}
2780:       MatAYPX(FDexp,-1.0,Bexp,SAME_NONZERO_PATTERN);
2781:       PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian\n");
2782:       if (flag) {MatView(FDexp,vstdout);}
2783:       if (vdraw) {              /* Always use contour for the difference */
2784:         PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);
2785:         MatView(FDexp,vdraw);
2786:         PetscViewerPopFormat(vdraw);
2787:       }
2788:       if (flag_contour) {PetscViewerPopFormat(vdraw);}
2789:       PetscViewerDestroy(&vdraw);
2790:       MatDestroy(&Bexp_mine);
2791:       MatDestroy(&FDexp);
2792:     }
2793:   }
2794:   {
2795:     PetscBool flag = PETSC_FALSE,flag_display = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_threshold = PETSC_FALSE;
2796:     PetscReal threshold_atol = PETSC_SQRT_MACHINE_EPSILON,threshold_rtol = 10*PETSC_SQRT_MACHINE_EPSILON;
2797:     PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_compare_coloring",NULL,NULL,&flag);
2798:     PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_compare_coloring_display",NULL,NULL,&flag_display);
2799:     PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_compare_coloring_draw",NULL,NULL,&flag_draw);
2800:     PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_compare_coloring_draw_contour",NULL,NULL,&flag_contour);
2801:     PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold",NULL,NULL,&flag_threshold);
2802:     if (flag_threshold) {
2803:       PetscOptionsGetReal(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_rtol",&threshold_rtol,NULL);
2804:       PetscOptionsGetReal(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_atol",&threshold_atol,NULL);
2805:     }
2806:     if (flag || flag_display || flag_draw || flag_contour || flag_threshold) {
2807:       Mat            Bfd;
2808:       PetscViewer    vdraw,vstdout;
2809:       MatColoring    coloring;
2810:       ISColoring     iscoloring;
2811:       MatFDColoring  matfdcoloring;
2812:       PetscErrorCode (*func)(SNES,Vec,Vec,void*);
2813:       void           *funcctx;
2814:       PetscReal      norm1,norm2,normmax;

2816:       MatDuplicate(B,MAT_DO_NOT_COPY_VALUES,&Bfd);
2817:       MatColoringCreate(Bfd,&coloring);
2818:       MatColoringSetType(coloring,MATCOLORINGSL);
2819:       MatColoringSetFromOptions(coloring);
2820:       MatColoringApply(coloring,&iscoloring);
2821:       MatColoringDestroy(&coloring);
2822:       MatFDColoringCreate(Bfd,iscoloring,&matfdcoloring);
2823:       MatFDColoringSetFromOptions(matfdcoloring);
2824:       MatFDColoringSetUp(Bfd,iscoloring,matfdcoloring);
2825:       ISColoringDestroy(&iscoloring);

2827:       /* This method of getting the function is currently unreliable since it doesn't work for DM local functions. */
2828:       SNESGetFunction(snes,NULL,&func,&funcctx);
2829:       MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))func,funcctx);
2830:       PetscObjectSetOptionsPrefix((PetscObject)matfdcoloring,((PetscObject)snes)->prefix);
2831:       PetscObjectAppendOptionsPrefix((PetscObject)matfdcoloring,"coloring_");
2832:       MatFDColoringSetFromOptions(matfdcoloring);
2833:       MatFDColoringApply(Bfd,matfdcoloring,X,snes);
2834:       MatFDColoringDestroy(&matfdcoloring);

2836:       PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&vstdout);
2837:       if (flag_draw || flag_contour) {
2838:         PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes),0,"Colored Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);
2839:         if (flag_contour) {PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);}
2840:       } else vdraw = NULL;
2841:       PetscViewerASCIIPrintf(vstdout,"Explicit preconditioning Jacobian\n");
2842:       if (flag_display) {MatView(B,vstdout);}
2843:       if (vdraw) {MatView(B,vdraw);}
2844:       PetscViewerASCIIPrintf(vstdout,"Colored Finite difference Jacobian\n");
2845:       if (flag_display) {MatView(Bfd,vstdout);}
2846:       if (vdraw) {MatView(Bfd,vdraw);}
2847:       MatAYPX(Bfd,-1.0,B,SAME_NONZERO_PATTERN);
2848:       MatNorm(Bfd,NORM_1,&norm1);
2849:       MatNorm(Bfd,NORM_FROBENIUS,&norm2);
2850:       MatNorm(Bfd,NORM_MAX,&normmax);
2851:       PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian, norm1=%g normFrob=%g normmax=%g\n",(double)norm1,(double)norm2,(double)normmax);
2852:       if (flag_display) {MatView(Bfd,vstdout);}
2853:       if (vdraw) {              /* Always use contour for the difference */
2854:         PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);
2855:         MatView(Bfd,vdraw);
2856:         PetscViewerPopFormat(vdraw);
2857:       }
2858:       if (flag_contour) {PetscViewerPopFormat(vdraw);}

2860:       if (flag_threshold) {
2861:         PetscInt bs,rstart,rend,i;
2862:         MatGetBlockSize(B,&bs);
2863:         MatGetOwnershipRange(B,&rstart,&rend);
2864:         for (i=rstart; i<rend; i++) {
2865:           const PetscScalar *ba,*ca;
2866:           const PetscInt    *bj,*cj;
2867:           PetscInt          bn,cn,j,maxentrycol = -1,maxdiffcol = -1,maxrdiffcol = -1;
2868:           PetscReal         maxentry = 0,maxdiff = 0,maxrdiff = 0;
2869:           MatGetRow(B,i,&bn,&bj,&ba);
2870:           MatGetRow(Bfd,i,&cn,&cj,&ca);
2871:           if (bn != cn) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_PLIB,"Unexpected different nonzero pattern in -snes_compare_coloring_threshold");
2872:           for (j=0; j<bn; j++) {
2873:             PetscReal rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j]));
2874:             if (PetscAbsScalar(ba[j]) > PetscAbs(maxentry)) {
2875:               maxentrycol = bj[j];
2876:               maxentry    = PetscRealPart(ba[j]);
2877:             }
2878:             if (PetscAbsScalar(ca[j]) > PetscAbs(maxdiff)) {
2879:               maxdiffcol = bj[j];
2880:               maxdiff    = PetscRealPart(ca[j]);
2881:             }
2882:             if (rdiff > maxrdiff) {
2883:               maxrdiffcol = bj[j];
2884:               maxrdiff    = rdiff;
2885:             }
2886:           }
2887:           if (maxrdiff > 1) {
2888:             PetscViewerASCIIPrintf(vstdout,"row %D (maxentry=%g at %D, maxdiff=%g at %D, maxrdiff=%g at %D):",i,(double)maxentry,maxentrycol,(double)maxdiff,maxdiffcol,(double)maxrdiff,maxrdiffcol);
2889:             for (j=0; j<bn; j++) {
2890:               PetscReal rdiff;
2891:               rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j]));
2892:               if (rdiff > 1) {
2893:                 PetscViewerASCIIPrintf(vstdout," (%D,%g:%g)",bj[j],(double)PetscRealPart(ba[j]),(double)PetscRealPart(ca[j]));
2894:               }
2895:             }
2896:             PetscViewerASCIIPrintf(vstdout,"\n",i,maxentry,maxdiff,maxrdiff);
2897:           }
2898:           MatRestoreRow(B,i,&bn,&bj,&ba);
2899:           MatRestoreRow(Bfd,i,&cn,&cj,&ca);
2900:         }
2901:       }
2902:       PetscViewerDestroy(&vdraw);
2903:       MatDestroy(&Bfd);
2904:     }
2905:   }
2906:   return(0);
2907: }

2909: /*MC
2910:     SNESJacobianFunction - Function used to convey the nonlinear Jacobian of the function to be solved by SNES

2912:      Synopsis:
2913:      #include "petscsnes.h"
2914:      PetscErrorCode SNESJacobianFunction(SNES snes,Vec x,Mat Amat,Mat Pmat,void *ctx);

2916: +  x - input vector
2917: .  Amat - the matrix that defines the (approximate) Jacobian
2918: .  Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.
2919: -  ctx - [optional] user-defined Jacobian context

2921:    Level: intermediate

2923: .seealso:   SNESSetFunction(), SNESGetFunction(), SNESSetJacobian(), SNESGetJacobian()
2924: M*/

2926: /*@C
2927:    SNESSetJacobian - Sets the function to compute Jacobian as well as the
2928:    location to store the matrix.

2930:    Logically Collective on SNES and Mat

2932:    Input Parameters:
2933: +  snes - the SNES context
2934: .  Amat - the matrix that defines the (approximate) Jacobian
2935: .  Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.
2936: .  J - Jacobian evaluation routine (if NULL then SNES retains any previously set value), see SNESJacobianFunction for details
2937: -  ctx - [optional] user-defined context for private data for the
2938:          Jacobian evaluation routine (may be NULL) (if NULL then SNES retains any previously set value)

2940:    Notes:
2941:    If the Amat matrix and Pmat matrix are different you must call MatAssemblyBegin/End() on
2942:    each matrix.

2944:    If you know the operator Amat has a null space you can use MatSetNullSpace() and MatSetTransposeNullSpace() to supply the null
2945:    space to Amat and the KSP solvers will automatically use that null space as needed during the solution process.

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

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

2953:    Level: beginner

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

2957: .seealso: KSPSetOperators(), SNESSetFunction(), MatMFFDComputeJacobian(), SNESComputeJacobianDefaultColor(), MatStructure, J,
2958:           SNESSetPicard(), SNESJacobianFunction
2959: @*/
2960: PetscErrorCode  SNESSetJacobian(SNES snes,Mat Amat,Mat Pmat,PetscErrorCode (*J)(SNES,Vec,Mat,Mat,void*),void *ctx)
2961: {
2963:   DM             dm;

2971:   SNESGetDM(snes,&dm);
2972:   DMSNESSetJacobian(dm,J,ctx);
2973:   if (Amat) {
2974:     PetscObjectReference((PetscObject)Amat);
2975:     MatDestroy(&snes->jacobian);

2977:     snes->jacobian = Amat;
2978:   }
2979:   if (Pmat) {
2980:     PetscObjectReference((PetscObject)Pmat);
2981:     MatDestroy(&snes->jacobian_pre);

2983:     snes->jacobian_pre = Pmat;
2984:   }
2985:   return(0);
2986: }

2988: /*@C
2989:    SNESGetJacobian - Returns the Jacobian matrix and optionally the user
2990:    provided context for evaluating the Jacobian.

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

2994:    Input Parameter:
2995: .  snes - the nonlinear solver context

2997:    Output Parameters:
2998: +  Amat - location to stash (approximate) Jacobian matrix (or NULL)
2999: .  Pmat - location to stash matrix used to compute the preconditioner (or NULL)
3000: .  J - location to put Jacobian function (or NULL), see SNESJacobianFunction for details on its calling sequence
3001: -  ctx - location to stash Jacobian ctx (or NULL)

3003:    Level: advanced

3005: .seealso: SNESSetJacobian(), SNESComputeJacobian(), SNESJacobianFunction, SNESGetFunction()
3006: @*/
3007: PetscErrorCode SNESGetJacobian(SNES snes,Mat *Amat,Mat *Pmat,PetscErrorCode (**J)(SNES,Vec,Mat,Mat,void*),void **ctx)
3008: {
3010:   DM             dm;
3011:   DMSNES         sdm;

3015:   if (Amat) *Amat = snes->jacobian;
3016:   if (Pmat) *Pmat = snes->jacobian_pre;
3017:   SNESGetDM(snes,&dm);
3018:   DMGetDMSNES(dm,&sdm);
3019:   if (J) *J = sdm->ops->computejacobian;
3020:   if (ctx) *ctx = sdm->jacobianctx;
3021:   return(0);
3022: }

3024: /*@
3025:    SNESSetUp - Sets up the internal data structures for the later use
3026:    of a nonlinear solver.

3028:    Collective on SNES

3030:    Input Parameters:
3031: .  snes - the SNES context

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

3040:    Level: advanced

3042: .keywords: SNES, nonlinear, setup

3044: .seealso: SNESCreate(), SNESSolve(), SNESDestroy()
3045: @*/
3046: PetscErrorCode  SNESSetUp(SNES snes)
3047: {
3049:   DM             dm;
3050:   DMSNES         sdm;
3051:   SNESLineSearch linesearch, pclinesearch;
3052:   void           *lsprectx,*lspostctx;
3053:   PetscErrorCode (*precheck)(SNESLineSearch,Vec,Vec,PetscBool*,void*);
3054:   PetscErrorCode (*postcheck)(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);
3055:   PetscErrorCode (*func)(SNES,Vec,Vec,void*);
3056:   Vec            f,fpc;
3057:   void           *funcctx;
3058:   PetscErrorCode (*jac)(SNES,Vec,Mat,Mat,void*);
3059:   void           *jacctx,*appctx;
3060:   Mat            j,jpre;

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

3066:   if (!((PetscObject)snes)->type_name) {
3067:     SNESSetType(snes,SNESNEWTONLS);
3068:   }

3070:   SNESGetFunction(snes,&snes->vec_func,NULL,NULL);

3072:   SNESGetDM(snes,&dm);
3073:   DMGetDMSNES(dm,&sdm);
3074:   if (!sdm->ops->computefunction) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Function never provided to SNES object");
3075:   if (!sdm->ops->computejacobian) {
3076:     DMSNESSetJacobian(dm,SNESComputeJacobianDefaultColor,NULL);
3077:   }
3078:   if (!snes->vec_func) {
3079:     DMCreateGlobalVector(dm,&snes->vec_func);
3080:   }

3082:   if (!snes->ksp) {
3083:     SNESGetKSP(snes, &snes->ksp);
3084:   }

3086:   if (!snes->linesearch) {
3087:     SNESGetLineSearch(snes, &snes->linesearch);
3088:   }
3089:   SNESLineSearchSetFunction(snes->linesearch,SNESComputeFunction);

3091:   if (snes->npc && (snes->npcside== PC_LEFT)) {
3092:     snes->mf          = PETSC_TRUE;
3093:     snes->mf_operator = PETSC_FALSE;
3094:   }

3096:   if (snes->npc) {
3097:     /* copy the DM over */
3098:     SNESGetDM(snes,&dm);
3099:     SNESSetDM(snes->npc,dm);

3101:     SNESGetFunction(snes,&f,&func,&funcctx);
3102:     VecDuplicate(f,&fpc);
3103:     SNESSetFunction(snes->npc,fpc,func,funcctx);
3104:     SNESGetJacobian(snes,&j,&jpre,&jac,&jacctx);
3105:     SNESSetJacobian(snes->npc,j,jpre,jac,jacctx);
3106:     SNESGetApplicationContext(snes,&appctx);
3107:     SNESSetApplicationContext(snes->npc,appctx);
3108:     VecDestroy(&fpc);

3110:     /* copy the function pointers over */
3111:     PetscObjectCopyFortranFunctionPointers((PetscObject)snes,(PetscObject)snes->npc);

3113:     /* default to 1 iteration */
3114:     SNESSetTolerances(snes->npc,0.0,0.0,0.0,1,snes->npc->max_funcs);
3115:     if (snes->npcside==PC_RIGHT) {
3116:       SNESSetNormSchedule(snes->npc,SNES_NORM_FINAL_ONLY);
3117:     } else {
3118:       SNESSetNormSchedule(snes->npc,SNES_NORM_NONE);
3119:     }
3120:     SNESSetFromOptions(snes->npc);

3122:     /* copy the line search context over */
3123:     SNESGetLineSearch(snes,&linesearch);
3124:     SNESGetLineSearch(snes->npc,&pclinesearch);
3125:     SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);
3126:     SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);
3127:     SNESLineSearchSetPreCheck(pclinesearch,precheck,lsprectx);
3128:     SNESLineSearchSetPostCheck(pclinesearch,postcheck,lspostctx);
3129:     PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)pclinesearch);
3130:   }
3131:   if (snes->mf) {
3132:     SNESSetUpMatrixFree_Private(snes, snes->mf_operator, snes->mf_version);
3133:   }
3134:   if (snes->ops->usercompute && !snes->user) {
3135:     (*snes->ops->usercompute)(snes,(void**)&snes->user);
3136:   }

3138:   snes->jac_iter = 0;
3139:   snes->pre_iter = 0;

3141:   if (snes->ops->setup) {
3142:     (*snes->ops->setup)(snes);
3143:   }

3145:   if (snes->npc && (snes->npcside== PC_LEFT)) {
3146:     if (snes->functype == SNES_FUNCTION_PRECONDITIONED) {
3147:       SNESGetLineSearch(snes,&linesearch);
3148:       SNESLineSearchSetFunction(linesearch,SNESComputeFunctionDefaultNPC);
3149:     }
3150:   }

3152:   snes->setupcalled = PETSC_TRUE;
3153:   return(0);
3154: }

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

3159:    Collective on SNES

3161:    Input Parameter:
3162: .  snes - iterative context obtained from SNESCreate()

3164:    Level: intermediate

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

3169: .keywords: SNES, destroy

3171: .seealso: SNESCreate(), SNESSetUp(), SNESSolve()
3172: @*/
3173: PetscErrorCode  SNESReset(SNES snes)
3174: {

3179:   if (snes->ops->userdestroy && snes->user) {
3180:     (*snes->ops->userdestroy)((void**)&snes->user);
3181:     snes->user = NULL;
3182:   }
3183:   if (snes->npc) {
3184:     SNESReset(snes->npc);
3185:   }

3187:   if (snes->ops->reset) {
3188:     (*snes->ops->reset)(snes);
3189:   }
3190:   if (snes->ksp) {
3191:     KSPReset(snes->ksp);
3192:   }

3194:   if (snes->linesearch) {
3195:     SNESLineSearchReset(snes->linesearch);
3196:   }

3198:   VecDestroy(&snes->vec_rhs);
3199:   VecDestroy(&snes->vec_sol);
3200:   VecDestroy(&snes->vec_sol_update);
3201:   VecDestroy(&snes->vec_func);
3202:   MatDestroy(&snes->jacobian);
3203:   MatDestroy(&snes->jacobian_pre);
3204:   VecDestroyVecs(snes->nwork,&snes->work);
3205:   VecDestroyVecs(snes->nvwork,&snes->vwork);

3207:   snes->alwayscomputesfinalresidual = PETSC_FALSE;

3209:   snes->nwork       = snes->nvwork = 0;
3210:   snes->setupcalled = PETSC_FALSE;
3211:   return(0);
3212: }

3214: /*@
3215:    SNESDestroy - Destroys the nonlinear solver context that was created
3216:    with SNESCreate().

3218:    Collective on SNES

3220:    Input Parameter:
3221: .  snes - the SNES context

3223:    Level: beginner

3225: .keywords: SNES, nonlinear, destroy

3227: .seealso: SNESCreate(), SNESSolve()
3228: @*/
3229: PetscErrorCode  SNESDestroy(SNES *snes)
3230: {

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

3238:   SNESReset((*snes));
3239:   SNESDestroy(&(*snes)->npc);

3241:   /* if memory was published with SAWs then destroy it */
3242:   PetscObjectSAWsViewOff((PetscObject)*snes);
3243:   if ((*snes)->ops->destroy) {(*((*snes))->ops->destroy)((*snes));}

3245:   if ((*snes)->dm) {DMCoarsenHookRemove((*snes)->dm,DMCoarsenHook_SNESVecSol,DMRestrictHook_SNESVecSol,*snes);}
3246:   DMDestroy(&(*snes)->dm);
3247:   KSPDestroy(&(*snes)->ksp);
3248:   SNESLineSearchDestroy(&(*snes)->linesearch);

3250:   PetscFree((*snes)->kspconvctx);
3251:   if ((*snes)->ops->convergeddestroy) {
3252:     (*(*snes)->ops->convergeddestroy)((*snes)->cnvP);
3253:   }
3254:   if ((*snes)->conv_malloc) {
3255:     PetscFree((*snes)->conv_hist);
3256:     PetscFree((*snes)->conv_hist_its);
3257:   }
3258:   SNESMonitorCancel((*snes));
3259:   PetscHeaderDestroy(snes);
3260:   return(0);
3261: }

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

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

3268:    Logically Collective on SNES

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

3275:    Options Database Keys:
3276: .    -snes_lag_preconditioner <lag>

3278:    Notes:
3279:    The default is 1
3280:    The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1
3281:    If  -1 is used before the very first nonlinear solve the preconditioner is still built because there is no previous preconditioner to use

3283:    Level: intermediate

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

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

3289: @*/
3290: PetscErrorCode  SNESSetLagPreconditioner(SNES snes,PetscInt lag)
3291: {
3294:   if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater");
3295:   if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0");
3297:   snes->lagpreconditioner = lag;
3298:   return(0);
3299: }

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

3304:    Logically Collective on SNES

3306:    Input Parameters:
3307: +  snes - the SNES context
3308: -  steps - the number of refinements to do, defaults to 0

3310:    Options Database Keys:
3311: .    -snes_grid_sequence <steps>

3313:    Level: intermediate

3315:    Notes:
3316:    Use SNESGetSolution() to extract the fine grid solution after grid sequencing.

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

3320: .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian(), SNESGetGridSequence()

3322: @*/
3323: PetscErrorCode  SNESSetGridSequence(SNES snes,PetscInt steps)
3324: {
3328:   snes->gridsequence = steps;
3329:   return(0);
3330: }

3332: /*@
3333:    SNESGetGridSequence - gets the number of steps of grid sequencing that SNES does

3335:    Logically Collective on SNES

3337:    Input Parameter:
3338: .  snes - the SNES context

3340:    Output Parameter:
3341: .  steps - the number of refinements to do, defaults to 0

3343:    Options Database Keys:
3344: .    -snes_grid_sequence <steps>

3346:    Level: intermediate

3348:    Notes:
3349:    Use SNESGetSolution() to extract the fine grid solution after grid sequencing.

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

3353: .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian(), SNESSetGridSequence()

3355: @*/
3356: PetscErrorCode  SNESGetGridSequence(SNES snes,PetscInt *steps)
3357: {
3360:   *steps = snes->gridsequence;
3361:   return(0);
3362: }

3364: /*@
3365:    SNESGetLagPreconditioner - Indicates how often the preconditioner is rebuilt

3367:    Not Collective

3369:    Input Parameter:
3370: .  snes - the SNES context

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

3376:    Options Database Keys:
3377: .    -snes_lag_preconditioner <lag>

3379:    Notes:
3380:    The default is 1
3381:    The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1

3383:    Level: intermediate

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

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

3389: @*/
3390: PetscErrorCode  SNESGetLagPreconditioner(SNES snes,PetscInt *lag)
3391: {
3394:   *lag = snes->lagpreconditioner;
3395:   return(0);
3396: }

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

3402:    Logically Collective on SNES

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

3409:    Options Database Keys:
3410: .    -snes_lag_jacobian <lag>

3412:    Notes:
3413:    The default is 1
3414:    The Jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1
3415:    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
3416:    at the next Newton step but never again (unless it is reset to another value)

3418:    Level: intermediate

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

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

3424: @*/
3425: PetscErrorCode  SNESSetLagJacobian(SNES snes,PetscInt lag)
3426: {
3429:   if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater");
3430:   if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0");
3432:   snes->lagjacobian = lag;
3433:   return(0);
3434: }

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

3439:    Not Collective

3441:    Input Parameter:
3442: .  snes - the SNES context

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

3448:    Options Database Keys:
3449: .    -snes_lag_jacobian <lag>

3451:    Notes:
3452:    The default is 1
3453:    The jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1

3455:    Level: intermediate

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

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

3461: @*/
3462: PetscErrorCode  SNESGetLagJacobian(SNES snes,PetscInt *lag)
3463: {
3466:   *lag = snes->lagjacobian;
3467:   return(0);
3468: }

3470: /*@
3471:    SNESSetLagJacobianPersists - Set whether or not the Jacobian lagging persists through multiple solves

3473:    Logically collective on SNES

3475:    Input Parameter:
3476: +  snes - the SNES context
3477: -   flg - jacobian lagging persists if true

3479:    Options Database Keys:
3480: .    -snes_lag_jacobian_persists <flg>

3482:    Notes:
3483:     This is useful both for nonlinear preconditioning, where it's appropriate to have the Jacobian be stale by
3484:    several solves, and for implicit time-stepping, where Jacobian lagging in the inner nonlinear solve over several
3485:    timesteps may present huge efficiency gains.

3487:    Level: developer

3489: .keywords: SNES, nonlinear, lag

3491: .seealso: SNESSetLagPreconditionerPersists(), SNESSetLagJacobian(), SNESGetLagJacobian(), SNESGetNPC()

3493: @*/
3494: PetscErrorCode  SNESSetLagJacobianPersists(SNES snes,PetscBool flg)
3495: {
3499:   snes->lagjac_persist = flg;
3500:   return(0);
3501: }

3503: /*@
3504:    SNESSetLagPreconditionerPersists - Set whether or not the preconditioner lagging persists through multiple solves

3506:    Logically Collective on SNES

3508:    Input Parameter:
3509: +  snes - the SNES context
3510: -   flg - preconditioner lagging persists if true

3512:    Options Database Keys:
3513: .    -snes_lag_jacobian_persists <flg>

3515:    Notes:
3516:     This is useful both for nonlinear preconditioning, where it's appropriate to have the preconditioner be stale
3517:    by several solves, and for implicit time-stepping, where preconditioner lagging in the inner nonlinear solve over
3518:    several timesteps may present huge efficiency gains.

3520:    Level: developer

3522: .keywords: SNES, nonlinear, lag

3524: .seealso: SNESSetLagJacobianPersists(), SNESSetLagJacobian(), SNESGetLagJacobian(), SNESGetNPC()

3526: @*/
3527: PetscErrorCode  SNESSetLagPreconditionerPersists(SNES snes,PetscBool flg)
3528: {
3532:   snes->lagpre_persist = flg;
3533:   return(0);
3534: }

3536: /*@
3537:    SNESSetForceIteration - force SNESSolve() to take at least one iteration regardless of the initial residual norm

3539:    Logically Collective on SNES

3541:    Input Parameters:
3542: +  snes - the SNES context
3543: -  force - PETSC_TRUE require at least one iteration

3545:    Options Database Keys:
3546: .    -snes_force_iteration <force> - Sets forcing an iteration

3548:    Notes:
3549:    This is used sometimes with TS to prevent TS from detecting a false steady state solution

3551:    Level: intermediate

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

3555: .seealso: SNESSetTrustRegionTolerance(), SNESSetDivergenceTolerance()
3556: @*/
3557: PetscErrorCode  SNESSetForceIteration(SNES snes,PetscBool force)
3558: {
3561:   snes->forceiteration = force;
3562:   return(0);
3563: }

3565: /*@
3566:    SNESGetForceIteration - Whether or not to force SNESSolve() take at least one iteration regardless of the initial residual norm

3568:    Logically Collective on SNES

3570:    Input Parameters:
3571: .  snes - the SNES context

3573:    Output Parameter:
3574: .  force - PETSC_TRUE requires at least one iteration.

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

3578:    Level: intermediate

3580: .seealso: SNESSetForceIteration(), SNESSetTrustRegionTolerance(), SNESSetDivergenceTolerance()
3581: @*/
3582: PetscErrorCode  SNESGetForceIteration(SNES snes,PetscBool *force)
3583: {
3586:   *force = snes->forceiteration;
3587:   return(0);
3588: }

3590: /*@
3591:    SNESSetTolerances - Sets various parameters used in convergence tests.

3593:    Logically Collective on SNES

3595:    Input Parameters:
3596: +  snes - the SNES context
3597: .  abstol - absolute convergence tolerance
3598: .  rtol - relative convergence tolerance
3599: .  stol -  convergence tolerance in terms of the norm of the change in the solution between steps,  || delta x || < stol*|| x ||
3600: .  maxit - maximum number of iterations
3601: -  maxf - maximum number of function evaluations (-1 indicates no limit)

3603:    Options Database Keys:
3604: +    -snes_atol <abstol> - Sets abstol
3605: .    -snes_rtol <rtol> - Sets rtol
3606: .    -snes_stol <stol> - Sets stol
3607: .    -snes_max_it <maxit> - Sets maxit
3608: -    -snes_max_funcs <maxf> - Sets maxf

3610:    Notes:
3611:    The default maximum number of iterations is 50.
3612:    The default maximum number of function evaluations is 1000.

3614:    Level: intermediate

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

3618: .seealso: SNESSetTrustRegionTolerance(), SNESSetDivergenceTolerance(), SNESSetForceIteration()
3619: @*/
3620: PetscErrorCode  SNESSetTolerances(SNES snes,PetscReal abstol,PetscReal rtol,PetscReal stol,PetscInt maxit,PetscInt maxf)
3621: {

3630:   if (abstol != PETSC_DEFAULT) {
3631:     if (abstol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %g must be non-negative",(double)abstol);
3632:     snes->abstol = abstol;
3633:   }
3634:   if (rtol != PETSC_DEFAULT) {
3635:     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",(double)rtol);
3636:     snes->rtol = rtol;
3637:   }
3638:   if (stol != PETSC_DEFAULT) {
3639:     if (stol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Step tolerance %g must be non-negative",(double)stol);
3640:     snes->stol = stol;
3641:   }
3642:   if (maxit != PETSC_DEFAULT) {
3643:     if (maxit < 0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxit);
3644:     snes->max_its = maxit;
3645:   }
3646:   if (maxf != PETSC_DEFAULT) {
3647:     if (maxf < -1) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of function evaluations %D must be -1 or nonnegative",maxf);
3648:     snes->max_funcs = maxf;
3649:   }
3650:   snes->tolerancesset = PETSC_TRUE;
3651:   return(0);
3652: }

3654: /*@
3655:    SNESSetDivergenceTolerance - Sets the divergence tolerance used for the SNES divergence test.

3657:    Logically Collective on SNES

3659:    Input Parameters:
3660: +  snes - the SNES context
3661: -  divtol - the divergence tolerance. Use -1 to deactivate the test.

3663:    Options Database Keys:
3664: +    -snes_divergence_tolerance <divtol> - Sets divtol

3666:    Notes:
3667:    The default divergence tolerance is 1e4.

3669:    Level: intermediate

3671: .keywords: SNES, nonlinear, set, divergence, tolerance

3673: .seealso: SNESSetTolerances(), SNESGetDivergenceTolerance
3674: @*/
3675: PetscErrorCode  SNESSetDivergenceTolerance(SNES snes,PetscReal divtol)
3676: {

3681:   if (divtol != PETSC_DEFAULT) {
3682:     snes->divtol = divtol;
3683:   }
3684:   else {
3685:     snes->divtol = 1.0e4;
3686:   }
3687:   return(0);
3688: }

3690: /*@
3691:    SNESGetTolerances - Gets various parameters used in convergence tests.

3693:    Not Collective

3695:    Input Parameters:
3696: +  snes - the SNES context
3697: .  atol - absolute convergence tolerance
3698: .  rtol - relative convergence tolerance
3699: .  stol -  convergence tolerance in terms of the norm
3700:            of the change in the solution between steps
3701: .  maxit - maximum number of iterations
3702: -  maxf - maximum number of function evaluations

3704:    Notes:
3705:    The user can specify NULL for any parameter that is not needed.

3707:    Level: intermediate

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

3711: .seealso: SNESSetTolerances()
3712: @*/
3713: PetscErrorCode  SNESGetTolerances(SNES snes,PetscReal *atol,PetscReal *rtol,PetscReal *stol,PetscInt *maxit,PetscInt *maxf)
3714: {
3717:   if (atol)  *atol  = snes->abstol;
3718:   if (rtol)  *rtol  = snes->rtol;
3719:   if (stol)  *stol  = snes->stol;
3720:   if (maxit) *maxit = snes->max_its;
3721:   if (maxf)  *maxf  = snes->max_funcs;
3722:   return(0);
3723: }

3725: /*@
3726:    SNESGetDivergenceTolerance - Gets divergence tolerance used in divergence test.

3728:    Not Collective

3730:    Input Parameters:
3731: +  snes - the SNES context
3732: -  divtol - divergence tolerance

3734:    Level: intermediate

3736: .keywords: SNES, nonlinear, get, divergence, tolerance

3738: .seealso: SNESSetDivergenceTolerance()
3739: @*/
3740: PetscErrorCode  SNESGetDivergenceTolerance(SNES snes,PetscReal *divtol)
3741: {
3744:   if (divtol) *divtol = snes->divtol;
3745:   return(0);
3746: }

3748: /*@
3749:    SNESSetTrustRegionTolerance - Sets the trust region parameter tolerance.

3751:    Logically Collective on SNES

3753:    Input Parameters:
3754: +  snes - the SNES context
3755: -  tol - tolerance

3757:    Options Database Key:
3758: .  -snes_trtol <tol> - Sets tol

3760:    Level: intermediate

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

3764: .seealso: SNESSetTolerances()
3765: @*/
3766: PetscErrorCode  SNESSetTrustRegionTolerance(SNES snes,PetscReal tol)
3767: {
3771:   snes->deltatol = tol;
3772:   return(0);
3773: }

3775: /*
3776:    Duplicate the lg monitors for SNES from KSP; for some reason with
3777:    dynamic libraries things don't work under Sun4 if we just use
3778:    macros instead of functions
3779: */
3780: PetscErrorCode  SNESMonitorLGResidualNorm(SNES snes,PetscInt it,PetscReal norm,void *ctx)
3781: {

3786:   KSPMonitorLGResidualNorm((KSP)snes,it,norm,ctx);
3787:   return(0);
3788: }

3790: PetscErrorCode  SNESMonitorLGCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *lgctx)
3791: {

3795:   KSPMonitorLGResidualNormCreate(comm,host,label,x,y,m,n,lgctx);
3796:   return(0);
3797: }

3799: PETSC_INTERN PetscErrorCode  SNESMonitorRange_Private(SNES,PetscInt,PetscReal*);

3801: PetscErrorCode  SNESMonitorLGRange(SNES snes,PetscInt n,PetscReal rnorm,void *monctx)
3802: {
3803:   PetscDrawLG      lg;
3804:   PetscErrorCode   ierr;
3805:   PetscReal        x,y,per;
3806:   PetscViewer      v = (PetscViewer)monctx;
3807:   static PetscReal prev; /* should be in the context */
3808:   PetscDraw        draw;

3812:   PetscViewerDrawGetDrawLG(v,0,&lg);
3813:   if (!n) {PetscDrawLGReset(lg);}
3814:   PetscDrawLGGetDraw(lg,&draw);
3815:   PetscDrawSetTitle(draw,"Residual norm");
3816:   x    = (PetscReal)n;
3817:   if (rnorm > 0.0) y = PetscLog10Real(rnorm);
3818:   else y = -15.0;
3819:   PetscDrawLGAddPoint(lg,&x,&y);
3820:   if (n < 20 || !(n % 5) || snes->reason) {
3821:     PetscDrawLGDraw(lg);
3822:     PetscDrawLGSave(lg);
3823:   }

3825:   PetscViewerDrawGetDrawLG(v,1,&lg);
3826:   if (!n) {PetscDrawLGReset(lg);}
3827:   PetscDrawLGGetDraw(lg,&draw);
3828:   PetscDrawSetTitle(draw,"% elemts > .2*max elemt");
3829:    SNESMonitorRange_Private(snes,n,&per);
3830:   x    = (PetscReal)n;
3831:   y    = 100.0*per;
3832:   PetscDrawLGAddPoint(lg,&x,&y);
3833:   if (n < 20 || !(n % 5) || snes->reason) {
3834:     PetscDrawLGDraw(lg);
3835:     PetscDrawLGSave(lg);
3836:   }

3838:   PetscViewerDrawGetDrawLG(v,2,&lg);
3839:   if (!n) {prev = rnorm;PetscDrawLGReset(lg);}
3840:   PetscDrawLGGetDraw(lg,&draw);
3841:   PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm");
3842:   x    = (PetscReal)n;
3843:   y    = (prev - rnorm)/prev;
3844:   PetscDrawLGAddPoint(lg,&x,&y);
3845:   if (n < 20 || !(n % 5) || snes->reason) {
3846:     PetscDrawLGDraw(lg);
3847:     PetscDrawLGSave(lg);
3848:   }

3850:   PetscViewerDrawGetDrawLG(v,3,&lg);
3851:   if (!n) {PetscDrawLGReset(lg);}
3852:   PetscDrawLGGetDraw(lg,&draw);
3853:   PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm*(% > .2 max)");
3854:   x    = (PetscReal)n;
3855:   y    = (prev - rnorm)/(prev*per);
3856:   if (n > 2) { /*skip initial crazy value */
3857:     PetscDrawLGAddPoint(lg,&x,&y);
3858:   }
3859:   if (n < 20 || !(n % 5) || snes->reason) {
3860:     PetscDrawLGDraw(lg);
3861:     PetscDrawLGSave(lg);
3862:   }
3863:   prev = rnorm;
3864:   return(0);
3865: }

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

3870:    Collective on SNES

3872:    Input Parameters:
3873: +  snes - nonlinear solver context obtained from SNESCreate()
3874: .  iter - iteration number
3875: -  rnorm - relative norm of the residual

3877:    Notes:
3878:    This routine is called by the SNES implementations.
3879:    It does not typically need to be called by the user.

3881:    Level: developer

3883: .seealso: SNESMonitorSet()
3884: @*/
3885: PetscErrorCode  SNESMonitor(SNES snes,PetscInt iter,PetscReal rnorm)
3886: {
3888:   PetscInt       i,n = snes->numbermonitors;

3891:   VecLockReadPush(snes->vec_sol);
3892:   for (i=0; i<n; i++) {
3893:     (*snes->monitor[i])(snes,iter,rnorm,snes->monitorcontext[i]);
3894:   }
3895:   VecLockReadPop(snes->vec_sol);
3896:   return(0);
3897: }

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

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

3904:      Synopsis:
3905:      #include <petscsnes.h>
3906: $    PetscErrorCode SNESMonitorFunction(SNES snes,PetscInt its, PetscReal norm,void *mctx)

3908: +    snes - the SNES context
3909: .    its - iteration number
3910: .    norm - 2-norm function value (may be estimated)
3911: -    mctx - [optional] monitoring context

3913:    Level: advanced

3915: .seealso:   SNESMonitorSet(), SNESMonitorGet()
3916: M*/

3918: /*@C
3919:    SNESMonitorSet - Sets an ADDITIONAL function that is to be used at every
3920:    iteration of the nonlinear solver to display the iteration's
3921:    progress.

3923:    Logically Collective on SNES

3925:    Input Parameters:
3926: +  snes - the SNES context
3927: .  f - the monitor function, see SNESMonitorFunction for the calling sequence
3928: .  mctx - [optional] user-defined context for private data for the
3929:           monitor routine (use NULL if no context is desired)
3930: -  monitordestroy - [optional] routine that frees monitor context
3931:           (may be NULL)

3933:    Options Database Keys:
3934: +    -snes_monitor        - sets SNESMonitorDefault()
3935: .    -snes_monitor_lg_residualnorm    - sets line graph monitor,
3936:                             uses SNESMonitorLGCreate()
3937: -    -snes_monitor_cancel - cancels all monitors that have
3938:                             been hardwired into a code by
3939:                             calls to SNESMonitorSet(), but
3940:                             does not cancel those set via
3941:                             the options database.

3943:    Notes:
3944:    Several different monitoring routines may be set by calling
3945:    SNESMonitorSet() multiple times; all will be called in the
3946:    order in which they were set.

3948:    Fortran Notes:
3949:     Only a single monitor function can be set for each SNES object

3951:    Level: intermediate

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

3955: .seealso: SNESMonitorDefault(), SNESMonitorCancel(), SNESMonitorFunction
3956: @*/
3957: PetscErrorCode  SNESMonitorSet(SNES snes,PetscErrorCode (*f)(SNES,PetscInt,PetscReal,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**))
3958: {
3959:   PetscInt       i;
3961:   PetscBool      identical;

3965:   for (i=0; i<snes->numbermonitors;i++) {
3966:     PetscMonitorCompare((PetscErrorCode (*)(void))f,mctx,monitordestroy,(PetscErrorCode (*)(void))snes->monitor[i],snes->monitorcontext[i],snes->monitordestroy[i],&identical);
3967:     if (identical) return(0);
3968:   }
3969:   if (snes->numbermonitors >= MAXSNESMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set");
3970:   snes->monitor[snes->numbermonitors]          = f;
3971:   snes->monitordestroy[snes->numbermonitors]   = monitordestroy;
3972:   snes->monitorcontext[snes->numbermonitors++] = (void*)mctx;
3973:   return(0);
3974: }

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

3979:    Logically Collective on SNES

3981:    Input Parameters:
3982: .  snes - the SNES context

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

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

3992:    Level: intermediate

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

3996: .seealso: SNESMonitorDefault(), SNESMonitorSet()
3997: @*/
3998: PetscErrorCode  SNESMonitorCancel(SNES snes)
3999: {
4001:   PetscInt       i;

4005:   for (i=0; i<snes->numbermonitors; i++) {
4006:     if (snes->monitordestroy[i]) {
4007:       (*snes->monitordestroy[i])(&snes->monitorcontext[i]);
4008:     }
4009:   }
4010:   snes->numbermonitors = 0;
4011:   return(0);
4012: }

4014: /*MC
4015:     SNESConvergenceTestFunction - functional form used for testing of convergence of nonlinear solver

4017:      Synopsis:
4018:      #include <petscsnes.h>
4019: $     PetscErrorCode SNESConvergenceTest(SNES snes,PetscInt it,PetscReal xnorm,PetscReal gnorm,PetscReal f,SNESConvergedReason *reason,void *cctx)

4021: +    snes - the SNES context
4022: .    it - current iteration (0 is the first and is before any Newton step)
4023: .    cctx - [optional] convergence context
4024: .    reason - reason for convergence/divergence
4025: .    xnorm - 2-norm of current iterate
4026: .    gnorm - 2-norm of current step
4027: -    f - 2-norm of function

4029:    Level: intermediate

4031: .seealso:   SNESSetConvergenceTest(), SNESGetConvergenceTest()
4032: M*/

4034: /*@C
4035:    SNESSetConvergenceTest - Sets the function that is to be used
4036:    to test for convergence of the nonlinear iterative solution.

4038:    Logically Collective on SNES

4040:    Input Parameters:
4041: +  snes - the SNES context
4042: .  SNESConvergenceTestFunction - routine to test for convergence
4043: .  cctx - [optional] context for private data for the convergence routine  (may be NULL)
4044: -  destroy - [optional] destructor for the context (may be NULL; PETSC_NULL_FUNCTION in Fortran)

4046:    Level: advanced

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

4050: .seealso: SNESConvergedDefault(), SNESConvergedSkip(), SNESConvergenceTestFunction
4051: @*/
4052: PetscErrorCode  SNESSetConvergenceTest(SNES snes,PetscErrorCode (*SNESConvergenceTestFunction)(SNES,PetscInt,PetscReal,PetscReal,PetscReal,SNESConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*))
4053: {

4058:   if (!SNESConvergenceTestFunction) SNESConvergenceTestFunction = SNESConvergedSkip;
4059:   if (snes->ops->convergeddestroy) {
4060:     (*snes->ops->convergeddestroy)(snes->cnvP);
4061:   }
4062:   snes->ops->converged        = SNESConvergenceTestFunction;
4063:   snes->ops->convergeddestroy = destroy;
4064:   snes->cnvP                  = cctx;
4065:   return(0);
4066: }

4068: /*@
4069:    SNESGetConvergedReason - Gets the reason the SNES iteration was stopped.

4071:    Not Collective

4073:    Input Parameter:
4074: .  snes - the SNES context

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

4080:    Options Database:
4081: .   -snes_converged_reason - prints the reason to standard out

4083:    Level: intermediate

4085:    Notes:
4086:     Should only be called after the call the SNESSolve() is complete, if it is called earlier it returns the value SNES__CONVERGED_ITERATING.

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

4090: .seealso: SNESSetConvergenceTest(), SNESSetConvergedReason(), SNESConvergedReason
4091: @*/
4092: PetscErrorCode SNESGetConvergedReason(SNES snes,SNESConvergedReason *reason)
4093: {
4097:   *reason = snes->reason;
4098:   return(0);
4099: }

4101: /*@
4102:    SNESSetConvergedReason - Sets the reason the SNES iteration was stopped.

4104:    Not Collective

4106:    Input Parameters:
4107: +  snes - the SNES context
4108: -  reason - negative value indicates diverged, positive value converged, see SNESConvergedReason or the
4109:             manual pages for the individual convergence tests for complete lists

4111:    Level: intermediate

4113: .keywords: SNES, nonlinear, set, convergence, test
4114: .seealso: SNESGetConvergedReason(), SNESSetConvergenceTest(), SNESConvergedReason
4115: @*/
4116: PetscErrorCode SNESSetConvergedReason(SNES snes,SNESConvergedReason reason)
4117: {
4120:   snes->reason = reason;
4121:   return(0);
4122: }

4124: /*@
4125:    SNESSetConvergenceHistory - Sets the array used to hold the convergence history.

4127:    Logically Collective on SNES

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

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

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

4145:    Level: intermediate

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

4149: .seealso: SNESGetConvergenceHistory()

4151: @*/
4152: PetscErrorCode  SNESSetConvergenceHistory(SNES snes,PetscReal a[],PetscInt its[],PetscInt na,PetscBool reset)
4153: {

4160:   if (!a) {
4161:     if (na == PETSC_DECIDE || na == PETSC_DEFAULT) na = 1000;
4162:     PetscCalloc1(na,&a);
4163:     PetscCalloc1(na,&its);

4165:     snes->conv_malloc = PETSC_TRUE;
4166:   }
4167:   snes->conv_hist       = a;
4168:   snes->conv_hist_its   = its;
4169:   snes->conv_hist_max   = na;
4170:   snes->conv_hist_len   = 0;
4171:   snes->conv_hist_reset = reset;
4172:   return(0);
4173: }

4175: #if defined(PETSC_HAVE_MATLAB_ENGINE)
4176: #include <engine.h>   /* MATLAB include file */
4177: #include <mex.h>      /* MATLAB include file */

4179: PETSC_EXTERN mxArray *SNESGetConvergenceHistoryMatlab(SNES snes)
4180: {
4181:   mxArray   *mat;
4182:   PetscInt  i;
4183:   PetscReal *ar;

4186:   mat = mxCreateDoubleMatrix(snes->conv_hist_len,1,mxREAL);
4187:   ar  = (PetscReal*) mxGetData(mat);
4188:   for (i=0; i<snes->conv_hist_len; i++) ar[i] = snes->conv_hist[i];
4189:   PetscFunctionReturn(mat);
4190: }
4191: #endif

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

4196:    Not Collective

4198:    Input Parameter:
4199: .  snes - iterative context obtained from SNESCreate()

4201:    Output Parameters:
4202: .  a   - array to hold history
4203: .  its - integer array holds the number of linear iterations (or
4204:          negative if not converged) for each solve.
4205: -  na  - size of a and its

4207:    Notes:
4208:     The calling sequence for this routine in Fortran is
4209: $   call SNESGetConvergenceHistory(SNES snes, integer na, integer ierr)

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

4215:    Level: intermediate

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

4219: .seealso: SNESSetConvergencHistory()

4221: @*/
4222: PetscErrorCode  SNESGetConvergenceHistory(SNES snes,PetscReal *a[],PetscInt *its[],PetscInt *na)
4223: {
4226:   if (a)   *a   = snes->conv_hist;
4227:   if (its) *its = snes->conv_hist_its;
4228:   if (na)  *na  = snes->conv_hist_len;
4229:   return(0);
4230: }

4232: /*@C
4233:   SNESSetUpdate - Sets the general-purpose update function called
4234:   at the beginning of every iteration of the nonlinear solve. Specifically
4235:   it is called just before the Jacobian is "evaluated".

4237:   Logically Collective on SNES

4239:   Input Parameters:
4240: . snes - The nonlinear solver context
4241: . func - The function

4243:   Calling sequence of func:
4244: . func (SNES snes, PetscInt step);

4246: . step - The current step of the iteration

4248:   Level: advanced

4250:   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()
4251:         This is not used by most users.

4253: .keywords: SNES, update

4255: .seealso SNESSetJacobian(), SNESSolve()
4256: @*/
4257: PetscErrorCode  SNESSetUpdate(SNES snes, PetscErrorCode (*func)(SNES, PetscInt))
4258: {
4261:   snes->ops->update = func;
4262:   return(0);
4263: }

4265: /*
4266:    SNESScaleStep_Private - Scales a step so that its length is less than the
4267:    positive parameter delta.

4269:     Input Parameters:
4270: +   snes - the SNES context
4271: .   y - approximate solution of linear system
4272: .   fnorm - 2-norm of current function
4273: -   delta - trust region size

4275:     Output Parameters:
4276: +   gpnorm - predicted function norm at the new point, assuming local
4277:     linearization.  The value is zero if the step lies within the trust
4278:     region, and exceeds zero otherwise.
4279: -   ynorm - 2-norm of the step

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

4285: .keywords: SNES, nonlinear, scale, step
4286: */
4287: PetscErrorCode SNESScaleStep_Private(SNES snes,Vec y,PetscReal *fnorm,PetscReal *delta,PetscReal *gpnorm,PetscReal *ynorm)
4288: {
4289:   PetscReal      nrm;
4290:   PetscScalar    cnorm;


4298:   VecNorm(y,NORM_2,&nrm);
4299:   if (nrm > *delta) {
4300:     nrm     = *delta/nrm;
4301:     *gpnorm = (1.0 - nrm)*(*fnorm);
4302:     cnorm   = nrm;
4303:     VecScale(y,cnorm);
4304:     *ynorm  = *delta;
4305:   } else {
4306:     *gpnorm = 0.0;
4307:     *ynorm  = nrm;
4308:   }
4309:   return(0);
4310: }

4312: /*@
4313:    SNESReasonView - Displays the reason a SNES solve converged or diverged to a viewer

4315:    Collective on SNES

4317:    Parameter:
4318: +  snes - iterative context obtained from SNESCreate()
4319: -  viewer - the viewer to display the reason


4322:    Options Database Keys:
4323: .  -snes_converged_reason - print reason for converged or diverged, also prints number of iterations

4325:    Level: beginner

4327: .keywords: SNES, solve, linear system

4329: .seealso: SNESCreate(), SNESSetUp(), SNESDestroy(), SNESSetTolerances(), SNESConvergedDefault()

4331: @*/
4332: PetscErrorCode  SNESReasonView(SNES snes,PetscViewer viewer)
4333: {
4334:   PetscViewerFormat format;
4335:   PetscBool         isAscii;
4336:   PetscErrorCode    ierr;

4339:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii);
4340:   if (isAscii) {
4341:     PetscViewerGetFormat(viewer, &format);
4342:     PetscViewerASCIIAddTab(viewer,((PetscObject)snes)->tablevel);
4343:     if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
4344:       DM                dm;
4345:       Vec               u;
4346:       PetscDS           prob;
4347:       PetscInt          Nf, f;
4348:       PetscErrorCode (**exactFuncs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
4349:       PetscReal         error;

4351:       SNESGetDM(snes, &dm);
4352:       SNESGetSolution(snes, &u);
4353:       DMGetDS(dm, &prob);
4354:       PetscDSGetNumFields(prob, &Nf);
4355:       PetscMalloc1(Nf, &exactFuncs);
4356:       for (f = 0; f < Nf; ++f) {PetscDSGetExactSolution(prob, f, &exactFuncs[f]);}
4357:       DMComputeL2Diff(dm, 0.0, exactFuncs, NULL, u, &error);
4358:       PetscFree(exactFuncs);
4359:       if (error < 1.0e-11) {PetscViewerASCIIPrintf(viewer, "L_2 Error: < 1.0e-11\n");}
4360:       else                 {PetscViewerASCIIPrintf(viewer, "L_2 Error: %g\n", error);}
4361:     }
4362:     if (snes->reason > 0) {
4363:       if (((PetscObject) snes)->prefix) {
4364:         PetscViewerASCIIPrintf(viewer,"Nonlinear %s solve converged due to %s iterations %D\n",((PetscObject) snes)->prefix,SNESConvergedReasons[snes->reason],snes->iter);
4365:       } else {
4366:         PetscViewerASCIIPrintf(viewer,"Nonlinear solve converged due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);
4367:       }
4368:     } else {
4369:       if (((PetscObject) snes)->prefix) {
4370:         PetscViewerASCIIPrintf(viewer,"Nonlinear %s solve did not converge due to %s iterations %D\n",((PetscObject) snes)->prefix,SNESConvergedReasons[snes->reason],snes->iter);
4371:       } else {
4372:         PetscViewerASCIIPrintf(viewer,"Nonlinear solve did not converge due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);
4373:       }
4374:     }
4375:     PetscViewerASCIISubtractTab(viewer,((PetscObject)snes)->tablevel);
4376:   }
4377:   return(0);
4378: }

4380: /*@C
4381:   SNESReasonViewFromOptions - Processes command line options to determine if/how a SNESReason is to be viewed.

4383:   Collective on SNES

4385:   Input Parameters:
4386: . snes   - the SNES object

4388:   Level: intermediate

4390: @*/
4391: PetscErrorCode SNESReasonViewFromOptions(SNES snes)
4392: {
4393:   PetscErrorCode    ierr;
4394:   PetscViewer       viewer;
4395:   PetscBool         flg;
4396:   static PetscBool  incall = PETSC_FALSE;
4397:   PetscViewerFormat format;

4400:   if (incall) return(0);
4401:   incall = PETSC_TRUE;
4402:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_converged_reason",&viewer,&format,&flg);
4403:   if (flg) {
4404:     PetscViewerPushFormat(viewer,format);
4405:     SNESReasonView(snes,viewer);
4406:     PetscViewerPopFormat(viewer);
4407:     PetscViewerDestroy(&viewer);
4408:   }
4409:   incall = PETSC_FALSE;
4410:   return(0);
4411: }

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

4417:    Collective on SNES

4419:    Input Parameters:
4420: +  snes - the SNES context
4421: .  b - the constant part of the equation F(x) = b, or NULL to use zero.
4422: -  x - the solution vector.

4424:    Notes:
4425:    The user should initialize the vector,x, with the initial guess
4426:    for the nonlinear solve prior to calling SNESSolve.  In particular,
4427:    to employ an initial guess of zero, the user should explicitly set
4428:    this vector to zero by calling VecSet().

4430:    Level: beginner

4432: .keywords: SNES, nonlinear, solve

4434: .seealso: SNESCreate(), SNESDestroy(), SNESSetFunction(), SNESSetJacobian(), SNESSetGridSequence(), SNESGetSolution()
4435: @*/
4436: PetscErrorCode  SNESSolve(SNES snes,Vec b,Vec x)
4437: {
4438:   PetscErrorCode    ierr;
4439:   PetscBool         flg;
4440:   PetscInt          grid;
4441:   Vec               xcreated = NULL;
4442:   DM                dm;


4451:   /* High level operations using the nonlinear solver */
4452:   {
4453:     PetscViewer       viewer;
4454:     PetscViewerFormat format;
4455:     PetscInt          num;
4456:     PetscBool         flg;
4457:     static PetscBool  incall = PETSC_FALSE;

4459:     if (!incall) {
4460:       /* Estimate the convergence rate of the discretization */
4461:       PetscOptionsGetViewer(PetscObjectComm((PetscObject) snes),((PetscObject)snes)->options, ((PetscObject) snes)->prefix, "-snes_convergence_estimate", &viewer, &format, &flg);
4462:       if (flg) {
4463:         PetscConvEst conv;
4464:         DM           dm;
4465:         PetscReal   *alpha; /* Convergence rate of the solution error for each field in the L_2 norm */
4466:         PetscInt     Nf;

4468:         incall = PETSC_TRUE;
4469:         SNESGetDM(snes, &dm);
4470:         DMGetNumFields(dm, &Nf);
4471:         PetscMalloc1(Nf, &alpha);
4472:         PetscConvEstCreate(PetscObjectComm((PetscObject) snes), &conv);
4473:         PetscConvEstSetSolver(conv, snes);
4474:         PetscConvEstSetFromOptions(conv);
4475:         PetscConvEstSetUp(conv);
4476:         PetscConvEstGetConvRate(conv, alpha);
4477:         PetscViewerPushFormat(viewer, format);
4478:         PetscConvEstRateView(conv, alpha, viewer);
4479:         PetscViewerPopFormat(viewer);
4480:         PetscViewerDestroy(&viewer);
4481:         PetscConvEstDestroy(&conv);
4482:         PetscFree(alpha);
4483:         incall = PETSC_FALSE;
4484:       }
4485:       /* Adaptively refine the initial grid */
4486:       num  = 1;
4487:       PetscOptionsGetInt(NULL, ((PetscObject) snes)->prefix, "-snes_adapt_initial", &num, &flg);
4488:       if (flg) {
4489:         DMAdaptor adaptor;

4491:         incall = PETSC_TRUE;
4492:         DMAdaptorCreate(PETSC_COMM_WORLD, &adaptor);
4493:         DMAdaptorSetSolver(adaptor, snes);
4494:         DMAdaptorSetSequenceLength(adaptor, num);
4495:         DMAdaptorSetFromOptions(adaptor);
4496:         DMAdaptorSetUp(adaptor);
4497:         DMAdaptorAdapt(adaptor, x, DM_ADAPTATION_INITIAL, &dm, &x);
4498:         DMAdaptorDestroy(&adaptor);
4499:         incall = PETSC_FALSE;
4500:       }
4501:       /* Use grid sequencing to adapt */
4502:       num  = 0;
4503:       PetscOptionsGetInt(NULL, ((PetscObject) snes)->prefix, "-snes_adapt_sequence", &num, NULL);
4504:       if (num) {
4505:         DMAdaptor adaptor;

4507:         incall = PETSC_TRUE;
4508:         DMAdaptorCreate(PETSC_COMM_WORLD, &adaptor);
4509:         DMAdaptorSetSolver(adaptor, snes);
4510:         DMAdaptorSetSequenceLength(adaptor, num);
4511:         DMAdaptorSetFromOptions(adaptor);
4512:         DMAdaptorSetUp(adaptor);
4513:         DMAdaptorAdapt(adaptor, x, DM_ADAPTATION_SEQUENTIAL, &dm, &x);
4514:         DMAdaptorDestroy(&adaptor);
4515:         incall = PETSC_FALSE;
4516:       }
4517:     }
4518:   }
4519:   if (!x) {
4520:     SNESGetDM(snes,&dm);
4521:     DMCreateGlobalVector(dm,&xcreated);
4522:     x    = xcreated;
4523:   }
4524:   SNESViewFromOptions(snes,NULL,"-snes_view_pre");

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

4529:     /* set solution vector */
4530:     if (!grid) {PetscObjectReference((PetscObject)x);}
4531:     VecDestroy(&snes->vec_sol);
4532:     snes->vec_sol = x;
4533:     SNESGetDM(snes,&dm);

4535:     /* set affine vector if provided */
4536:     if (b) { PetscObjectReference((PetscObject)b); }
4537:     VecDestroy(&snes->vec_rhs);
4538:     snes->vec_rhs = b;

4540:     if (snes->vec_rhs && (snes->vec_func == snes->vec_rhs)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Right hand side vector cannot be function vector");
4541:     if (snes->vec_func == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be function vector");
4542:     if (snes->vec_rhs  == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be right hand side vector");
4543:     if (!snes->vec_sol_update /* && snes->vec_sol */) {
4544:       VecDuplicate(snes->vec_sol,&snes->vec_sol_update);
4545:       PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->vec_sol_update);
4546:     }
4547:     DMShellSetGlobalVector(dm,snes->vec_sol);
4548:     SNESSetUp(snes);

4550:     if (!grid) {
4551:       if (snes->ops->computeinitialguess) {
4552:         (*snes->ops->computeinitialguess)(snes,snes->vec_sol,snes->initialguessP);
4553:       }
4554:     }

4556:     if (snes->conv_hist_reset) snes->conv_hist_len = 0;
4557:     if (snes->counters_reset) {snes->nfuncs = 0; snes->linear_its = 0; snes->numFailures = 0;}

4559:     PetscLogEventBegin(SNES_Solve,snes,0,0,0);
4560:     (*snes->ops->solve)(snes);
4561:     PetscLogEventEnd(SNES_Solve,snes,0,0,0);
4562:     if (!snes->reason) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
4563:     snes->domainerror = PETSC_FALSE; /* clear the flag if it has been set */

4565:     if (snes->lagjac_persist) snes->jac_iter += snes->iter;
4566:     if (snes->lagpre_persist) snes->pre_iter += snes->iter;

4568:     PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_test_local_min",NULL,NULL,&flg);
4569:     if (flg && !PetscPreLoadingOn) { SNESTestLocalMin(snes); }
4570:     SNESReasonViewFromOptions(snes);

4572:     if (snes->errorifnotconverged && snes->reason < 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged");
4573:     if (snes->reason < 0) break;
4574:     if (grid <  snes->gridsequence) {
4575:       DM  fine;
4576:       Vec xnew;
4577:       Mat interp;

4579:       DMRefine(snes->dm,PetscObjectComm((PetscObject)snes),&fine);
4580:       if (!fine) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"DMRefine() did not perform any refinement, cannot continue grid sequencing");
4581:       DMCreateInterpolation(snes->dm,fine,&interp,NULL);
4582:       DMCreateGlobalVector(fine,&xnew);
4583:       MatInterpolate(interp,x,xnew);
4584:       DMInterpolate(snes->dm,interp,fine);
4585:       MatDestroy(&interp);
4586:       x    = xnew;

4588:       SNESReset(snes);
4589:       SNESSetDM(snes,fine);
4590:       SNESResetFromOptions(snes);
4591:       DMDestroy(&fine);
4592:       PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));
4593:     }
4594:   }
4595:   SNESViewFromOptions(snes,NULL,"-snes_view");
4596:   VecViewFromOptions(snes->vec_sol,(PetscObject)snes,"-snes_view_solution");

4598:   VecDestroy(&xcreated);
4599:   PetscObjectSAWsBlock((PetscObject)snes);
4600:   return(0);
4601: }

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

4605: /*@C
4606:    SNESSetType - Sets the method for the nonlinear solver.

4608:    Collective on SNES

4610:    Input Parameters:
4611: +  snes - the SNES context
4612: -  type - a known method

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

4618:    Notes:
4619:    See "petsc/include/petscsnes.h" for available methods (for instance)
4620: +    SNESNEWTONLS - Newton's method with line search
4621:      (systems of nonlinear equations)
4622: .    SNESNEWTONTR - Newton's method with trust region
4623:      (systems of nonlinear equations)

4625:   Normally, it is best to use the SNESSetFromOptions() command and then
4626:   set the SNES solver type from the options database rather than by using
4627:   this routine.  Using the options database provides the user with
4628:   maximum flexibility in evaluating the many nonlinear solvers.
4629:   The SNESSetType() routine is provided for those situations where it
4630:   is necessary to set the nonlinear solver independently of the command
4631:   line or options database.  This might be the case, for example, when
4632:   the choice of solver changes during the execution of the program,
4633:   and the user's application is taking responsibility for choosing the
4634:   appropriate method.

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

4640:   Level: intermediate

4642: .keywords: SNES, set, type

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

4646: @*/
4647: PetscErrorCode  SNESSetType(SNES snes,SNESType type)
4648: {
4649:   PetscErrorCode ierr,(*r)(SNES);
4650:   PetscBool      match;


4656:   PetscObjectTypeCompare((PetscObject)snes,type,&match);
4657:   if (match) return(0);

4659:    PetscFunctionListFind(SNESList,type,&r);
4660:   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested SNES type %s",type);
4661:   /* Destroy the previous private SNES context */
4662:   if (snes->ops->destroy) {
4663:     (*(snes)->ops->destroy)(snes);
4664:     snes->ops->destroy = NULL;
4665:   }
4666:   /* Reinitialize function pointers in SNESOps structure */
4667:   snes->ops->setup          = 0;
4668:   snes->ops->solve          = 0;
4669:   snes->ops->view           = 0;
4670:   snes->ops->setfromoptions = 0;
4671:   snes->ops->destroy        = 0;
4672:   /* Call the SNESCreate_XXX routine for this particular Nonlinear solver */
4673:   snes->setupcalled = PETSC_FALSE;

4675:   PetscObjectChangeTypeName((PetscObject)snes,type);
4676:   (*r)(snes);
4677:   return(0);
4678: }

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

4683:    Not Collective

4685:    Input Parameter:
4686: .  snes - nonlinear solver context

4688:    Output Parameter:
4689: .  type - SNES method (a character string)

4691:    Level: intermediate

4693: .keywords: SNES, nonlinear, get, type, name
4694: @*/
4695: PetscErrorCode  SNESGetType(SNES snes,SNESType *type)
4696: {
4700:   *type = ((PetscObject)snes)->type_name;
4701:   return(0);
4702: }

4704: /*@
4705:   SNESSetSolution - Sets the solution vector for use by the SNES routines.

4707:   Logically Collective on SNES and Vec

4709:   Input Parameters:
4710: + snes - the SNES context obtained from SNESCreate()
4711: - u    - the solution vector

4713:   Level: beginner

4715: .keywords: SNES, set, solution
4716: @*/
4717: PetscErrorCode SNESSetSolution(SNES snes, Vec u)
4718: {
4719:   DM             dm;

4725:   PetscObjectReference((PetscObject) u);
4726:   VecDestroy(&snes->vec_sol);

4728:   snes->vec_sol = u;

4730:   SNESGetDM(snes, &dm);
4731:   DMShellSetGlobalVector(dm, u);
4732:   return(0);
4733: }

4735: /*@
4736:    SNESGetSolution - Returns the vector where the approximate solution is
4737:    stored. This is the fine grid solution when using SNESSetGridSequence().

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

4741:    Input Parameter:
4742: .  snes - the SNES context

4744:    Output Parameter:
4745: .  x - the solution

4747:    Level: intermediate

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

4751: .seealso:  SNESGetSolutionUpdate(), SNESGetFunction()
4752: @*/
4753: PetscErrorCode  SNESGetSolution(SNES snes,Vec *x)
4754: {
4758:   *x = snes->vec_sol;
4759:   return(0);
4760: }

4762: /*@
4763:    SNESGetSolutionUpdate - Returns the vector where the solution update is
4764:    stored.

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

4768:    Input Parameter:
4769: .  snes - the SNES context

4771:    Output Parameter:
4772: .  x - the solution update

4774:    Level: advanced

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

4778: .seealso: SNESGetSolution(), SNESGetFunction()
4779: @*/
4780: PetscErrorCode  SNESGetSolutionUpdate(SNES snes,Vec *x)
4781: {
4785:   *x = snes->vec_sol_update;
4786:   return(0);
4787: }

4789: /*@C
4790:    SNESGetFunction - Returns the vector where the function is stored.

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

4794:    Input Parameter:
4795: .  snes - the SNES context

4797:    Output Parameter:
4798: +  r - the vector that is used to store residuals (or NULL if you don't want it)
4799: .  f - the function (or NULL if you don't want it); see SNESFunction for calling sequence details
4800: -  ctx - the function context (or NULL if you don't want it)

4802:    Level: advanced

4804:     Notes: The vector r DOES NOT, in general contain the current value of the SNES nonlinear function

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

4808: .seealso: SNESSetFunction(), SNESGetSolution(), SNESFunction
4809: @*/
4810: PetscErrorCode  SNESGetFunction(SNES snes,Vec *r,PetscErrorCode (**f)(SNES,Vec,Vec,void*),void **ctx)
4811: {
4813:   DM             dm;

4817:   if (r) {
4818:     if (!snes->vec_func) {
4819:       if (snes->vec_rhs) {
4820:         VecDuplicate(snes->vec_rhs,&snes->vec_func);
4821:       } else if (snes->vec_sol) {
4822:         VecDuplicate(snes->vec_sol,&snes->vec_func);
4823:       } else if (snes->dm) {
4824:         DMCreateGlobalVector(snes->dm,&snes->vec_func);
4825:       }
4826:     }
4827:     *r = snes->vec_func;
4828:   }
4829:   SNESGetDM(snes,&dm);
4830:   DMSNESGetFunction(dm,f,ctx);
4831:   return(0);
4832: }

4834: /*@C
4835:    SNESGetNGS - Returns the NGS function and context.

4837:    Input Parameter:
4838: .  snes - the SNES context

4840:    Output Parameter:
4841: +  f - the function (or NULL) see SNESNGSFunction for details
4842: -  ctx    - the function context (or NULL)

4844:    Level: advanced

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

4848: .seealso: SNESSetNGS(), SNESGetFunction()
4849: @*/

4851: PetscErrorCode SNESGetNGS (SNES snes, PetscErrorCode (**f)(SNES, Vec, Vec, void*), void ** ctx)
4852: {
4854:   DM             dm;

4858:   SNESGetDM(snes,&dm);
4859:   DMSNESGetNGS(dm,f,ctx);
4860:   return(0);
4861: }

4863: /*@C
4864:    SNESSetOptionsPrefix - Sets the prefix used for searching for all
4865:    SNES options in the database.

4867:    Logically Collective on SNES

4869:    Input Parameter:
4870: +  snes - the SNES context
4871: -  prefix - the prefix to prepend to all option names

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

4877:    Level: advanced

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

4881: .seealso: SNESSetFromOptions()
4882: @*/
4883: PetscErrorCode  SNESSetOptionsPrefix(SNES snes,const char prefix[])
4884: {

4889:   PetscObjectSetOptionsPrefix((PetscObject)snes,prefix);
4890:   if (!snes->ksp) {SNESGetKSP(snes,&snes->ksp);}
4891:   if (snes->linesearch) {
4892:     SNESGetLineSearch(snes,&snes->linesearch);
4893:     PetscObjectSetOptionsPrefix((PetscObject)snes->linesearch,prefix);
4894:   }
4895:   KSPSetOptionsPrefix(snes->ksp,prefix);
4896:   return(0);
4897: }

4899: /*@C
4900:    SNESAppendOptionsPrefix - Appends to the prefix used for searching for all
4901:    SNES options in the database.

4903:    Logically Collective on SNES

4905:    Input Parameters:
4906: +  snes - the SNES context
4907: -  prefix - the prefix to prepend to all option names

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

4913:    Level: advanced

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

4917: .seealso: SNESGetOptionsPrefix()
4918: @*/
4919: PetscErrorCode  SNESAppendOptionsPrefix(SNES snes,const char prefix[])
4920: {

4925:   PetscObjectAppendOptionsPrefix((PetscObject)snes,prefix);
4926:   if (!snes->ksp) {SNESGetKSP(snes,&snes->ksp);}
4927:   if (snes->linesearch) {
4928:     SNESGetLineSearch(snes,&snes->linesearch);
4929:     PetscObjectAppendOptionsPrefix((PetscObject)snes->linesearch,prefix);
4930:   }
4931:   KSPAppendOptionsPrefix(snes->ksp,prefix);
4932:   return(0);
4933: }

4935: /*@C
4936:    SNESGetOptionsPrefix - Sets the prefix used for searching for all
4937:    SNES options in the database.

4939:    Not Collective

4941:    Input Parameter:
4942: .  snes - the SNES context

4944:    Output Parameter:
4945: .  prefix - pointer to the prefix string used

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

4951:    Level: advanced

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

4955: .seealso: SNESAppendOptionsPrefix()
4956: @*/
4957: PetscErrorCode  SNESGetOptionsPrefix(SNES snes,const char *prefix[])
4958: {

4963:   PetscObjectGetOptionsPrefix((PetscObject)snes,prefix);
4964:   return(0);
4965: }


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

4971:    Not collective

4973:    Input Parameters:
4974: +  name_solver - name of a new user-defined solver
4975: -  routine_create - routine to create method context

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

4980:    Sample usage:
4981: .vb
4982:    SNESRegister("my_solver",MySolverCreate);
4983: .ve

4985:    Then, your solver can be chosen with the procedural interface via
4986: $     SNESSetType(snes,"my_solver")
4987:    or at runtime via the option
4988: $     -snes_type my_solver

4990:    Level: advanced

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

4994: .keywords: SNES, nonlinear, register

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

4998:   Level: advanced
4999: @*/
5000: PetscErrorCode  SNESRegister(const char sname[],PetscErrorCode (*function)(SNES))
5001: {

5005:   SNESInitializePackage();
5006:   PetscFunctionListAdd(&SNESList,sname,function);
5007:   return(0);
5008: }

5010: PetscErrorCode  SNESTestLocalMin(SNES snes)
5011: {
5013:   PetscInt       N,i,j;
5014:   Vec            u,uh,fh;
5015:   PetscScalar    value;
5016:   PetscReal      norm;

5019:   SNESGetSolution(snes,&u);
5020:   VecDuplicate(u,&uh);
5021:   VecDuplicate(u,&fh);

5023:   /* currently only works for sequential */
5024:   PetscPrintf(PETSC_COMM_WORLD,"Testing FormFunction() for local min\n");
5025:   VecGetSize(u,&N);
5026:   for (i=0; i<N; i++) {
5027:     VecCopy(u,uh);
5028:     PetscPrintf(PETSC_COMM_WORLD,"i = %D\n",i);
5029:     for (j=-10; j<11; j++) {
5030:       value = PetscSign(j)*PetscExpReal(PetscAbs(j)-10.0);
5031:       VecSetValue(uh,i,value,ADD_VALUES);
5032:       SNESComputeFunction(snes,uh,fh);
5033:       VecNorm(fh,NORM_2,&norm);
5034:       PetscPrintf(PETSC_COMM_WORLD,"       j norm %D %18.16e\n",j,norm);
5035:       value = -value;
5036:       VecSetValue(uh,i,value,ADD_VALUES);
5037:     }
5038:   }
5039:   VecDestroy(&uh);
5040:   VecDestroy(&fh);
5041:   return(0);
5042: }

5044: /*@
5045:    SNESKSPSetUseEW - Sets SNES use Eisenstat-Walker method for
5046:    computing relative tolerance for linear solvers within an inexact
5047:    Newton method.

5049:    Logically Collective on SNES

5051:    Input Parameters:
5052: +  snes - SNES context
5053: -  flag - PETSC_TRUE or PETSC_FALSE

5055:     Options Database:
5056: +  -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence
5057: .  -snes_ksp_ew_version ver - version of  Eisenstat-Walker method
5058: .  -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0
5059: .  -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax
5060: .  -snes_ksp_ew_gamma <gamma> - Sets gamma
5061: .  -snes_ksp_ew_alpha <alpha> - Sets alpha
5062: .  -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2
5063: -  -snes_ksp_ew_threshold <threshold> - Sets threshold

5065:    Notes:
5066:    Currently, the default is to use a constant relative tolerance for
5067:    the inner linear solvers.  Alternatively, one can use the
5068:    Eisenstat-Walker method, where the relative convergence tolerance
5069:    is reset at each Newton iteration according progress of the nonlinear
5070:    solver.

5072:    Level: advanced

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

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

5080: .seealso: SNESKSPGetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW()
5081: @*/
5082: PetscErrorCode  SNESKSPSetUseEW(SNES snes,PetscBool flag)
5083: {
5087:   snes->ksp_ewconv = flag;
5088:   return(0);
5089: }

5091: /*@
5092:    SNESKSPGetUseEW - Gets if SNES is using Eisenstat-Walker method
5093:    for computing relative tolerance for linear solvers within an
5094:    inexact Newton method.

5096:    Not Collective

5098:    Input Parameter:
5099: .  snes - SNES context

5101:    Output Parameter:
5102: .  flag - PETSC_TRUE or PETSC_FALSE

5104:    Notes:
5105:    Currently, the default is to use a constant relative tolerance for
5106:    the inner linear solvers.  Alternatively, one can use the
5107:    Eisenstat-Walker method, where the relative convergence tolerance
5108:    is reset at each Newton iteration according progress of the nonlinear
5109:    solver.

5111:    Level: advanced

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

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

5119: .seealso: SNESKSPSetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW()
5120: @*/
5121: PetscErrorCode  SNESKSPGetUseEW(SNES snes, PetscBool  *flag)
5122: {
5126:   *flag = snes->ksp_ewconv;
5127:   return(0);
5128: }

5130: /*@
5131:    SNESKSPSetParametersEW - Sets parameters for Eisenstat-Walker
5132:    convergence criteria for the linear solvers within an inexact
5133:    Newton method.

5135:    Logically Collective on SNES

5137:    Input Parameters:
5138: +    snes - SNES context
5139: .    version - version 1, 2 (default is 2) or 3
5140: .    rtol_0 - initial relative tolerance (0 <= rtol_0 < 1)
5141: .    rtol_max - maximum relative tolerance (0 <= rtol_max < 1)
5142: .    gamma - multiplicative factor for version 2 rtol computation
5143:              (0 <= gamma2 <= 1)
5144: .    alpha - power for version 2 rtol computation (1 < alpha <= 2)
5145: .    alpha2 - power for safeguard
5146: -    threshold - threshold for imposing safeguard (0 < threshold < 1)

5148:    Note:
5149:    Version 3 was contributed by Luis Chacon, June 2006.

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

5153:    Level: advanced

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

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

5162: .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPGetParametersEW()
5163: @*/
5164: PetscErrorCode  SNESKSPSetParametersEW(SNES snes,PetscInt version,PetscReal rtol_0,PetscReal rtol_max,PetscReal gamma,PetscReal alpha,PetscReal alpha2,PetscReal threshold)
5165: {
5166:   SNESKSPEW *kctx;

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

5180:   if (version != PETSC_DEFAULT)   kctx->version   = version;
5181:   if (rtol_0 != PETSC_DEFAULT)    kctx->rtol_0    = rtol_0;
5182:   if (rtol_max != PETSC_DEFAULT)  kctx->rtol_max  = rtol_max;
5183:   if (gamma != PETSC_DEFAULT)     kctx->gamma     = gamma;
5184:   if (alpha != PETSC_DEFAULT)     kctx->alpha     = alpha;
5185:   if (alpha2 != PETSC_DEFAULT)    kctx->alpha2    = alpha2;
5186:   if (threshold != PETSC_DEFAULT) kctx->threshold = threshold;

5188:   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);
5189:   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",(double)kctx->rtol_0);
5190:   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",(double)kctx->rtol_max);
5191:   if (kctx->gamma < 0.0 || kctx->gamma > 1.0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= gamma (%g) <= 1.0\n",(double)kctx->gamma);
5192:   if (kctx->alpha <= 1.0 || kctx->alpha > 2.0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"1.0 < alpha (%g) <= 2.0\n",(double)kctx->alpha);
5193:   if (kctx->threshold <= 0.0 || kctx->threshold >= 1.0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 < threshold (%g) < 1.0\n",(double)kctx->threshold);
5194:   return(0);
5195: }

5197: /*@
5198:    SNESKSPGetParametersEW - Gets parameters for Eisenstat-Walker
5199:    convergence criteria for the linear solvers within an inexact
5200:    Newton method.

5202:    Not Collective

5204:    Input Parameters:
5205:      snes - SNES context

5207:    Output Parameters:
5208: +    version - version 1, 2 (default is 2) or 3
5209: .    rtol_0 - initial relative tolerance (0 <= rtol_0 < 1)
5210: .    rtol_max - maximum relative tolerance (0 <= rtol_max < 1)
5211: .    gamma - multiplicative factor for version 2 rtol computation (0 <= gamma2 <= 1)
5212: .    alpha - power for version 2 rtol computation (1 < alpha <= 2)
5213: .    alpha2 - power for safeguard
5214: -    threshold - threshold for imposing safeguard (0 < threshold < 1)

5216:    Level: advanced

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

5220: .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPSetParametersEW()
5221: @*/
5222: PetscErrorCode  SNESKSPGetParametersEW(SNES snes,PetscInt *version,PetscReal *rtol_0,PetscReal *rtol_max,PetscReal *gamma,PetscReal *alpha,PetscReal *alpha2,PetscReal *threshold)
5223: {
5224:   SNESKSPEW *kctx;

5228:   kctx = (SNESKSPEW*)snes->kspconvctx;
5229:   if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing");
5230:   if (version)   *version   = kctx->version;
5231:   if (rtol_0)    *rtol_0    = kctx->rtol_0;
5232:   if (rtol_max)  *rtol_max  = kctx->rtol_max;
5233:   if (gamma)     *gamma     = kctx->gamma;
5234:   if (alpha)     *alpha     = kctx->alpha;
5235:   if (alpha2)    *alpha2    = kctx->alpha2;
5236:   if (threshold) *threshold = kctx->threshold;
5237:   return(0);
5238: }

5240:  PetscErrorCode KSPPreSolve_SNESEW(KSP ksp, Vec b, Vec x, SNES snes)
5241: {
5243:   SNESKSPEW      *kctx = (SNESKSPEW*)snes->kspconvctx;
5244:   PetscReal      rtol  = PETSC_DEFAULT,stol;

5247:   if (!snes->ksp_ewconv) return(0);
5248:   if (!snes->iter) {
5249:     rtol = kctx->rtol_0; /* first time in, so use the original user rtol */
5250:     VecNorm(snes->vec_func,NORM_2,&kctx->norm_first);
5251:   }
5252:   else {
5253:     if (kctx->version == 1) {
5254:       rtol = (snes->norm - kctx->lresid_last)/kctx->norm_last;
5255:       if (rtol < 0.0) rtol = -rtol;
5256:       stol = PetscPowReal(kctx->rtol_last,kctx->alpha2);
5257:       if (stol > kctx->threshold) rtol = PetscMax(rtol,stol);
5258:     } else if (kctx->version == 2) {
5259:       rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha);
5260:       stol = kctx->gamma * PetscPowReal(kctx->rtol_last,kctx->alpha);
5261:       if (stol > kctx->threshold) rtol = PetscMax(rtol,stol);
5262:     } else if (kctx->version == 3) { /* contributed by Luis Chacon, June 2006. */
5263:       rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha);
5264:       /* safeguard: avoid sharp decrease of rtol */
5265:       stol = kctx->gamma*PetscPowReal(kctx->rtol_last,kctx->alpha);
5266:       stol = PetscMax(rtol,stol);
5267:       rtol = PetscMin(kctx->rtol_0,stol);
5268:       /* safeguard: avoid oversolving */
5269:       stol = kctx->gamma*(kctx->norm_first*snes->rtol)/snes->norm;
5270:       stol = PetscMax(rtol,stol);
5271:       rtol = PetscMin(kctx->rtol_0,stol);
5272:     } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 or 3 are supported: %D",kctx->version);
5273:   }
5274:   /* safeguard: avoid rtol greater than one */
5275:   rtol = PetscMin(rtol,kctx->rtol_max);
5276:   KSPSetTolerances(ksp,rtol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
5277:   PetscInfo3(snes,"iter %D, Eisenstat-Walker (version %D) KSP rtol=%g\n",snes->iter,kctx->version,(double)rtol);
5278:   return(0);
5279: }

5281: PetscErrorCode KSPPostSolve_SNESEW(KSP ksp, Vec b, Vec x, SNES snes)
5282: {
5284:   SNESKSPEW      *kctx = (SNESKSPEW*)snes->kspconvctx;
5285:   PCSide         pcside;
5286:   Vec            lres;

5289:   if (!snes->ksp_ewconv) return(0);
5290:   KSPGetTolerances(ksp,&kctx->rtol_last,0,0,0);
5291:   kctx->norm_last = snes->norm;
5292:   if (kctx->version == 1) {
5293:     PC        pc;
5294:     PetscBool isNone;

5296:     KSPGetPC(ksp, &pc);
5297:     PetscObjectTypeCompare((PetscObject) pc, PCNONE, &isNone);
5298:     KSPGetPCSide(ksp,&pcside);
5299:      if (pcside == PC_RIGHT || isNone) { /* XXX Should we also test KSP_UNPRECONDITIONED_NORM ? */
5300:       /* KSP residual is true linear residual */
5301:       KSPGetResidualNorm(ksp,&kctx->lresid_last);
5302:     } else {
5303:       /* KSP residual is preconditioned residual */
5304:       /* compute true linear residual norm */
5305:       VecDuplicate(b,&lres);
5306:       MatMult(snes->jacobian,x,lres);
5307:       VecAYPX(lres,-1.0,b);
5308:       VecNorm(lres,NORM_2,&kctx->lresid_last);
5309:       VecDestroy(&lres);
5310:     }
5311:   }
5312:   return(0);
5313: }

5315: /*@
5316:    SNESGetKSP - Returns the KSP context for a SNES solver.

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

5320:    Input Parameter:
5321: .  snes - the SNES context

5323:    Output Parameter:
5324: .  ksp - the KSP context

5326:    Notes:
5327:    The user can then directly manipulate the KSP context to set various
5328:    options, etc.  Likewise, the user can then extract and manipulate the
5329:    PC contexts as well.

5331:    Level: beginner

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

5335: .seealso: KSPGetPC(), SNESCreate(), KSPCreate(), SNESSetKSP()
5336: @*/
5337: PetscErrorCode  SNESGetKSP(SNES snes,KSP *ksp)
5338: {


5345:   if (!snes->ksp) {
5346:     PetscBool monitor = PETSC_FALSE;

5348:     KSPCreate(PetscObjectComm((PetscObject)snes),&snes->ksp);
5349:     PetscObjectIncrementTabLevel((PetscObject)snes->ksp,(PetscObject)snes,1);
5350:     PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->ksp);

5352:     KSPSetPreSolve(snes->ksp,(PetscErrorCode (*)(KSP,Vec,Vec,void*))KSPPreSolve_SNESEW,snes);
5353:     KSPSetPostSolve(snes->ksp,(PetscErrorCode (*)(KSP,Vec,Vec,void*))KSPPostSolve_SNESEW,snes);

5355:     PetscOptionsGetBool(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-ksp_monitor_snes",&monitor,NULL);
5356:     if (monitor) {
5357:       KSPMonitorSet(snes->ksp,KSPMonitorSNES,snes,NULL);
5358:     }
5359:     monitor = PETSC_FALSE;
5360:     PetscOptionsGetBool(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-ksp_monitor_snes_lg",&monitor,NULL);
5361:     if (monitor) {
5362:       PetscObject *objs;
5363:       KSPMonitorSNESLGResidualNormCreate(PetscObjectComm((PetscObject)snes),NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,600,600,&objs);
5364:       objs[0] = (PetscObject) snes;
5365:       KSPMonitorSet(snes->ksp,(PetscErrorCode (*)(KSP,PetscInt,PetscReal,void*))KSPMonitorSNESLGResidualNorm,objs,(PetscErrorCode (*)(void**))KSPMonitorSNESLGResidualNormDestroy);
5366:     }
5367:     PetscObjectSetOptions((PetscObject)snes->ksp,((PetscObject)snes)->options);
5368:   }
5369:   *ksp = snes->ksp;
5370:   return(0);
5371: }


5374:  #include <petsc/private/dmimpl.h>
5375: /*@
5376:    SNESSetDM - Sets the DM that may be used by some nonlinear solvers or their underlying preconditioners

5378:    Logically Collective on SNES

5380:    Input Parameters:
5381: +  snes - the nonlinear solver context
5382: -  dm - the dm, cannot be NULL

5384:    Notes:
5385:    A DM can only be used for solving one problem at a time because information about the problem is stored on the DM,
5386:    even when not using interfaces like DMSNESSetFunction().  Use DMClone() to get a distinct DM when solving different
5387:    problems using the same function space.

5389:    Level: intermediate

5391: .seealso: SNESGetDM(), KSPSetDM(), KSPGetDM()
5392: @*/
5393: PetscErrorCode  SNESSetDM(SNES snes,DM dm)
5394: {
5396:   KSP            ksp;
5397:   DMSNES         sdm;

5402:   PetscObjectReference((PetscObject)dm);
5403:   if (snes->dm) {               /* Move the DMSNES context over to the new DM unless the new DM already has one */
5404:     if (snes->dm->dmsnes && !dm->dmsnes) {
5405:       DMCopyDMSNES(snes->dm,dm);
5406:       DMGetDMSNES(snes->dm,&sdm);
5407:       if (sdm->originaldm == snes->dm) sdm->originaldm = dm; /* Grant write privileges to the replacement DM */
5408:     }
5409:     DMCoarsenHookRemove(snes->dm,DMCoarsenHook_SNESVecSol,DMRestrictHook_SNESVecSol,snes);
5410:     DMDestroy(&snes->dm);
5411:   }
5412:   snes->dm     = dm;
5413:   snes->dmAuto = PETSC_FALSE;

5415:   SNESGetKSP(snes,&ksp);
5416:   KSPSetDM(ksp,dm);
5417:   KSPSetDMActive(ksp,PETSC_FALSE);
5418:   if (snes->npc) {
5419:     SNESSetDM(snes->npc, snes->dm);
5420:     SNESSetNPCSide(snes,snes->npcside);
5421:   }
5422:   return(0);
5423: }

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

5428:    Not Collective but DM obtained is parallel on SNES

5430:    Input Parameter:
5431: . snes - the preconditioner context

5433:    Output Parameter:
5434: .  dm - the dm

5436:    Level: intermediate

5438: .seealso: SNESSetDM(), KSPSetDM(), KSPGetDM()
5439: @*/
5440: PetscErrorCode  SNESGetDM(SNES snes,DM *dm)
5441: {

5446:   if (!snes->dm) {
5447:     DMShellCreate(PetscObjectComm((PetscObject)snes),&snes->dm);
5448:     snes->dmAuto = PETSC_TRUE;
5449:   }
5450:   *dm = snes->dm;
5451:   return(0);
5452: }

5454: /*@
5455:   SNESSetNPC - Sets the nonlinear preconditioner to be used.

5457:   Collective on SNES

5459:   Input Parameters:
5460: + snes - iterative context obtained from SNESCreate()
5461: - pc   - the preconditioner object

5463:   Notes:
5464:   Use SNESGetNPC() to retrieve the preconditioner context (for example,
5465:   to configure it using the API).

5467:   Level: developer

5469: .keywords: SNES, set, precondition
5470: .seealso: SNESGetNPC(), SNESHasNPC()
5471: @*/
5472: PetscErrorCode SNESSetNPC(SNES snes, SNES pc)
5473: {

5480:   PetscObjectReference((PetscObject) pc);
5481:   SNESDestroy(&snes->npc);
5482:   snes->npc = pc;
5483:   PetscLogObjectParent((PetscObject)snes, (PetscObject)snes->npc);
5484:   return(0);
5485: }

5487: /*@
5488:   SNESGetNPC - Creates a nonlinear preconditioning solver (SNES) to be used to precondition the nonlinear solver.

5490:   Not Collective; but any changes to the obtained SNES object must be applied collectively

5492:   Input Parameter:
5493: . snes - iterative context obtained from SNESCreate()

5495:   Output Parameter:
5496: . pc - preconditioner context

5498:   Notes:
5499:     If a SNES was previously set with SNESSetNPC() then that SNES is returned.

5501:     The (preconditioner) SNES returned automatically inherits the same nonlinear function and Jacobian supplied to the original
5502:     SNES during SNESSetUp()

5504:   Level: developer

5506: .keywords: SNES, get, preconditioner
5507: .seealso: SNESSetNPC(), SNESHasNPC(), SNES, SNESCreate()
5508: @*/
5509: PetscErrorCode SNESGetNPC(SNES snes, SNES *pc)
5510: {
5512:   const char     *optionsprefix;

5517:   if (!snes->npc) {
5518:     SNESCreate(PetscObjectComm((PetscObject)snes),&snes->npc);
5519:     PetscObjectIncrementTabLevel((PetscObject)snes->npc,(PetscObject)snes,1);
5520:     PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->npc);
5521:     SNESGetOptionsPrefix(snes,&optionsprefix);
5522:     SNESSetOptionsPrefix(snes->npc,optionsprefix);
5523:     SNESAppendOptionsPrefix(snes->npc,"npc_");
5524:     SNESSetCountersReset(snes->npc,PETSC_FALSE);
5525:   }
5526:   *pc = snes->npc;
5527:   return(0);
5528: }

5530: /*@
5531:   SNESHasNPC - Returns whether a nonlinear preconditioner exists

5533:   Not Collective

5535:   Input Parameter:
5536: . snes - iterative context obtained from SNESCreate()

5538:   Output Parameter:
5539: . has_npc - whether the SNES has an NPC or not

5541:   Level: developer

5543: .keywords: SNES, has, preconditioner
5544: .seealso: SNESSetNPC(), SNESGetNPC()
5545: @*/
5546: PetscErrorCode SNESHasNPC(SNES snes, PetscBool *has_npc)
5547: {
5550:   *has_npc = (PetscBool) (snes->npc ? PETSC_TRUE : PETSC_FALSE);
5551:   return(0);
5552: }

5554: /*@
5555:     SNESSetNPCSide - Sets the preconditioning side.

5557:     Logically Collective on SNES

5559:     Input Parameter:
5560: .   snes - iterative context obtained from SNESCreate()

5562:     Output Parameter:
5563: .   side - the preconditioning side, where side is one of
5564: .vb
5565:       PC_LEFT - left preconditioning
5566:       PC_RIGHT - right preconditioning (default for most nonlinear solvers)
5567: .ve

5569:     Options Database Keys:
5570: .   -snes_pc_side <right,left>

5572:     Notes:
5573:     SNESNRICHARDSON and SNESNCG only support left preconditioning.

5575:     Level: intermediate

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

5579: .seealso: SNESGetNPCSide(), KSPSetPCSide()
5580: @*/
5581: PetscErrorCode  SNESSetNPCSide(SNES snes,PCSide side)
5582: {
5586:   snes->npcside= side;
5587:   return(0);
5588: }

5590: /*@
5591:     SNESGetNPCSide - Gets the preconditioning side.

5593:     Not Collective

5595:     Input Parameter:
5596: .   snes - iterative context obtained from SNESCreate()

5598:     Output Parameter:
5599: .   side - the preconditioning side, where side is one of
5600: .vb
5601:       PC_LEFT - left preconditioning
5602:       PC_RIGHT - right preconditioning (default for most nonlinear solvers)
5603: .ve

5605:     Level: intermediate

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

5609: .seealso: SNESSetNPCSide(), KSPGetPCSide()
5610: @*/
5611: PetscErrorCode  SNESGetNPCSide(SNES snes,PCSide *side)
5612: {
5616:   *side = snes->npcside;
5617:   return(0);
5618: }

5620: /*@
5621:   SNESSetLineSearch - Sets the linesearch on the SNES instance.

5623:   Collective on SNES

5625:   Input Parameters:
5626: + snes - iterative context obtained from SNESCreate()
5627: - linesearch   - the linesearch object

5629:   Notes:
5630:   Use SNESGetLineSearch() to retrieve the preconditioner context (for example,
5631:   to configure it using the API).

5633:   Level: developer

5635: .keywords: SNES, set, linesearch
5636: .seealso: SNESGetLineSearch()
5637: @*/
5638: PetscErrorCode SNESSetLineSearch(SNES snes, SNESLineSearch linesearch)
5639: {

5646:   PetscObjectReference((PetscObject) linesearch);
5647:   SNESLineSearchDestroy(&snes->linesearch);

5649:   snes->linesearch = linesearch;

5651:   PetscLogObjectParent((PetscObject)snes, (PetscObject)snes->linesearch);
5652:   return(0);
5653: }

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

5659:   Not Collective

5661:   Input Parameter:
5662: . snes - iterative context obtained from SNESCreate()

5664:   Output Parameter:
5665: . linesearch - linesearch context

5667:   Level: beginner

5669: .keywords: SNES, get, linesearch
5670: .seealso: SNESSetLineSearch(), SNESLineSearchCreate()
5671: @*/
5672: PetscErrorCode SNESGetLineSearch(SNES snes, SNESLineSearch *linesearch)
5673: {
5675:   const char     *optionsprefix;

5680:   if (!snes->linesearch) {
5681:     SNESGetOptionsPrefix(snes, &optionsprefix);
5682:     SNESLineSearchCreate(PetscObjectComm((PetscObject)snes), &snes->linesearch);
5683:     SNESLineSearchSetSNES(snes->linesearch, snes);
5684:     SNESLineSearchAppendOptionsPrefix(snes->linesearch, optionsprefix);
5685:     PetscObjectIncrementTabLevel((PetscObject) snes->linesearch, (PetscObject) snes, 1);
5686:     PetscLogObjectParent((PetscObject)snes, (PetscObject)snes->linesearch);
5687:   }
5688:   *linesearch = snes->linesearch;
5689:   return(0);
5690: }

5692: #if defined(PETSC_HAVE_MATLAB_ENGINE)
5693: #include <mex.h>

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

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

5700:    Collective on SNES

5702:    Input Parameters:
5703: +  snes - the SNES context
5704: -  x - input vector

5706:    Output Parameter:
5707: .  y - function vector, as set by SNESSetFunction()

5709:    Notes:
5710:    SNESComputeFunction() is typically used within nonlinear solvers
5711:    implementations, so most users would not generally call this routine
5712:    themselves.

5714:    Level: developer

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

5718: .seealso: SNESSetFunction(), SNESGetFunction()
5719: */
5720: PetscErrorCode  SNESComputeFunction_Matlab(SNES snes,Vec x,Vec y, void *ctx)
5721: {
5722:   PetscErrorCode    ierr;
5723:   SNESMatlabContext *sctx = (SNESMatlabContext*)ctx;
5724:   int               nlhs  = 1,nrhs = 5;
5725:   mxArray           *plhs[1],*prhs[5];
5726:   long long int     lx = 0,ly = 0,ls = 0;


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

5737:   PetscMemcpy(&ls,&snes,sizeof(snes));
5738:   PetscMemcpy(&lx,&x,sizeof(x));
5739:   PetscMemcpy(&ly,&y,sizeof(x));
5740:   prhs[0] = mxCreateDoubleScalar((double)ls);
5741:   prhs[1] = mxCreateDoubleScalar((double)lx);
5742:   prhs[2] = mxCreateDoubleScalar((double)ly);
5743:   prhs[3] = mxCreateString(sctx->funcname);
5744:   prhs[4] = sctx->ctx;
5745:   mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeFunctionInternal");
5746:   mxGetScalar(plhs[0]);
5747:   mxDestroyArray(prhs[0]);
5748:   mxDestroyArray(prhs[1]);
5749:   mxDestroyArray(prhs[2]);
5750:   mxDestroyArray(prhs[3]);
5751:   mxDestroyArray(plhs[0]);
5752:   return(0);
5753: }

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

5760:    Logically Collective on SNES

5762:    Input Parameters:
5763: +  snes - the SNES context
5764: .  r - vector to store function value
5765: -  f - function evaluation routine

5767:    Notes:
5768:    The Newton-like methods typically solve linear systems of the form
5769: $      f'(x) x = -f(x),
5770:    where f'(x) denotes the Jacobian matrix and f(x) is the function.

5772:    Level: beginner

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

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

5778: .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction()
5779: */
5780: PetscErrorCode  SNESSetFunctionMatlab(SNES snes,Vec r,const char *f,mxArray *ctx)
5781: {
5782:   PetscErrorCode    ierr;
5783:   SNESMatlabContext *sctx;

5786:   /* currently sctx is memory bleed */
5787:   PetscNew(&sctx);
5788:   PetscStrallocpy(f,&sctx->funcname);
5789:   /*
5790:      This should work, but it doesn't
5791:   sctx->ctx = ctx;
5792:   mexMakeArrayPersistent(sctx->ctx);
5793:   */
5794:   sctx->ctx = mxDuplicateArray(ctx);
5795:   SNESSetFunction(snes,r,SNESComputeFunction_Matlab,sctx);
5796:   return(0);
5797: }

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

5802:    Collective on SNES

5804:    Input Parameters:
5805: +  snes - the SNES context
5806: .  x - input vector
5807: .  A, B - the matrices
5808: -  ctx - user context

5810:    Level: developer

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

5814: .seealso: SNESSetFunction(), SNESGetFunction()
5815: @*/
5816: PetscErrorCode  SNESComputeJacobian_Matlab(SNES snes,Vec x,Mat A,Mat B,void *ctx)
5817: {
5818:   PetscErrorCode    ierr;
5819:   SNESMatlabContext *sctx = (SNESMatlabContext*)ctx;
5820:   int               nlhs  = 2,nrhs = 6;
5821:   mxArray           *plhs[2],*prhs[6];
5822:   long long int     lx = 0,lA = 0,ls = 0, lB = 0;


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

5830:   PetscMemcpy(&ls,&snes,sizeof(snes));
5831:   PetscMemcpy(&lx,&x,sizeof(x));
5832:   PetscMemcpy(&lA,A,sizeof(x));
5833:   PetscMemcpy(&lB,B,sizeof(x));
5834:   prhs[0] = mxCreateDoubleScalar((double)ls);
5835:   prhs[1] = mxCreateDoubleScalar((double)lx);
5836:   prhs[2] = mxCreateDoubleScalar((double)lA);
5837:   prhs[3] = mxCreateDoubleScalar((double)lB);
5838:   prhs[4] = mxCreateString(sctx->funcname);
5839:   prhs[5] = sctx->ctx;
5840:   mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeJacobianInternal");
5841:   mxGetScalar(plhs[0]);
5842:   mxDestroyArray(prhs[0]);
5843:   mxDestroyArray(prhs[1]);
5844:   mxDestroyArray(prhs[2]);
5845:   mxDestroyArray(prhs[3]);
5846:   mxDestroyArray(prhs[4]);
5847:   mxDestroyArray(plhs[0]);
5848:   mxDestroyArray(plhs[1]);
5849:   return(0);
5850: }

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

5857:    Logically Collective on SNES

5859:    Input Parameters:
5860: +  snes - the SNES context
5861: .  A,B - Jacobian matrices
5862: .  J - function evaluation routine
5863: -  ctx - user context

5865:    Level: developer

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

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

5871: .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction(), J
5872: */
5873: PetscErrorCode  SNESSetJacobianMatlab(SNES snes,Mat A,Mat B,const char *J,mxArray *ctx)
5874: {
5875:   PetscErrorCode    ierr;
5876:   SNESMatlabContext *sctx;

5879:   /* currently sctx is memory bleed */
5880:   PetscNew(&sctx);
5881:   PetscStrallocpy(J,&sctx->funcname);
5882:   /*
5883:      This should work, but it doesn't
5884:   sctx->ctx = ctx;
5885:   mexMakeArrayPersistent(sctx->ctx);
5886:   */
5887:   sctx->ctx = mxDuplicateArray(ctx);
5888:   SNESSetJacobian(snes,A,B,SNESComputeJacobian_Matlab,sctx);
5889:   return(0);
5890: }

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

5895:    Collective on SNES

5897: .seealso: SNESSetFunction(), SNESGetFunction()
5898: @*/
5899: PetscErrorCode  SNESMonitor_Matlab(SNES snes,PetscInt it, PetscReal fnorm, void *ctx)
5900: {
5901:   PetscErrorCode    ierr;
5902:   SNESMatlabContext *sctx = (SNESMatlabContext*)ctx;
5903:   int               nlhs  = 1,nrhs = 6;
5904:   mxArray           *plhs[1],*prhs[6];
5905:   long long int     lx = 0,ls = 0;
5906:   Vec               x  = snes->vec_sol;


5911:   PetscMemcpy(&ls,&snes,sizeof(snes));
5912:   PetscMemcpy(&lx,&x,sizeof(x));
5913:   prhs[0] = mxCreateDoubleScalar((double)ls);
5914:   prhs[1] = mxCreateDoubleScalar((double)it);
5915:   prhs[2] = mxCreateDoubleScalar((double)fnorm);
5916:   prhs[3] = mxCreateDoubleScalar((double)lx);
5917:   prhs[4] = mxCreateString(sctx->funcname);
5918:   prhs[5] = sctx->ctx;
5919:   mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESMonitorInternal");
5920:   mxGetScalar(plhs[0]);
5921:   mxDestroyArray(prhs[0]);
5922:   mxDestroyArray(prhs[1]);
5923:   mxDestroyArray(prhs[2]);
5924:   mxDestroyArray(prhs[3]);
5925:   mxDestroyArray(prhs[4]);
5926:   mxDestroyArray(plhs[0]);
5927:   return(0);
5928: }

5930: /*
5931:    SNESMonitorSetMatlab - Sets the monitor function from MATLAB

5933:    Level: developer

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

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

5939: .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction()
5940: */
5941: PetscErrorCode  SNESMonitorSetMatlab(SNES snes,const char *f,mxArray *ctx)
5942: {
5943:   PetscErrorCode    ierr;
5944:   SNESMatlabContext *sctx;

5947:   /* currently sctx is memory bleed */
5948:   PetscNew(&sctx);
5949:   PetscStrallocpy(f,&sctx->funcname);
5950:   /*
5951:      This should work, but it doesn't
5952:   sctx->ctx = ctx;
5953:   mexMakeArrayPersistent(sctx->ctx);
5954:   */
5955:   sctx->ctx = mxDuplicateArray(ctx);
5956:   SNESMonitorSet(snes,SNESMonitor_Matlab,sctx,NULL);
5957:   return(0);
5958: }

5960: #endif