Actual source code: snes.c
2: #include <private/snesimpl.h> /*I "petscsnes.h" I*/
4: PetscBool SNESRegisterAllCalled = PETSC_FALSE;
5: PetscFList SNESList = PETSC_NULL;
7: /* Logging support */
8: PetscClassId SNES_CLASSID;
9: PetscLogEvent SNES_Solve, SNES_LineSearch, SNES_FunctionEval, SNES_JacobianEval, SNES_GSEval;
13: /*
14: Translates from a SNES call to a DM call in computing a Jacobian
15: */
16: PetscErrorCode SNESDMComputeJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
17: {
19: DM dm;
22: SNESGetDM(snes,&dm);
23: DMComputeJacobian(dm,X,*J,*B,flag);
24: return(0);
25: }
29: /*@
30: SNESSetErrorIfNotConverged - Causes SNESSolve() to generate an error if the solver has not converged.
32: Logically Collective on SNES
34: Input Parameters:
35: + snes - iterative context obtained from SNESCreate()
36: - flg - PETSC_TRUE indicates you want the error generated
38: Options database keys:
39: . -snes_error_if_not_converged : this takes an optional truth value (0/1/no/yes/true/false)
41: Level: intermediate
43: Notes:
44: Normally PETSc continues if a linear solver fails to converge, you can call SNESGetConvergedReason() after a SNESSolve()
45: to determine if it has converged.
47: .keywords: SNES, set, initial guess, nonzero
49: .seealso: SNESGetErrorIfNotConverged(), KSPGetErrorIfNotConverged(), KSPSetErrorIFNotConverged()
50: @*/
51: PetscErrorCode SNESSetErrorIfNotConverged(SNES snes,PetscBool flg)
52: {
56: snes->errorifnotconverged = flg;
58: return(0);
59: }
63: /*@
64: SNESGetErrorIfNotConverged - Will SNESSolve() generate an error if the solver does not converge?
66: Not Collective
68: Input Parameter:
69: . snes - iterative context obtained from SNESCreate()
71: Output Parameter:
72: . flag - PETSC_TRUE if it will generate an error, else PETSC_FALSE
74: Level: intermediate
76: .keywords: SNES, set, initial guess, nonzero
78: .seealso: SNESSetErrorIfNotConverged(), KSPGetErrorIfNotConverged(), KSPSetErrorIFNotConverged()
79: @*/
80: PetscErrorCode SNESGetErrorIfNotConverged(SNES snes,PetscBool *flag)
81: {
85: *flag = snes->errorifnotconverged;
86: return(0);
87: }
91: /*@
92: SNESSetFunctionDomainError - tells SNES that the input vector to your FormFunction is not
93: in the functions domain. For example, negative pressure.
95: Logically Collective on SNES
97: Input Parameters:
98: . SNES - the SNES context
100: Level: advanced
102: .keywords: SNES, view
104: .seealso: SNESCreate(), SNESSetFunction()
105: @*/
106: PetscErrorCode SNESSetFunctionDomainError(SNES snes)
107: {
110: snes->domainerror = PETSC_TRUE;
111: return(0);
112: }
116: /*@C
117: SNESView - Prints the SNES data structure.
119: Collective on SNES
121: Input Parameters:
122: + SNES - the SNES context
123: - viewer - visualization context
125: Options Database Key:
126: . -snes_view - Calls SNESView() at end of SNESSolve()
128: Notes:
129: The available visualization contexts include
130: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
131: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
132: output where only the first processor opens
133: the file. All other processors send their
134: data to the first processor to print.
136: The user can open an alternative visualization context with
137: PetscViewerASCIIOpen() - output to a specified file.
139: Level: beginner
141: .keywords: SNES, view
143: .seealso: PetscViewerASCIIOpen()
144: @*/
145: PetscErrorCode SNESView(SNES snes,PetscViewer viewer)
146: {
147: SNESKSPEW *kctx;
148: PetscErrorCode ierr;
149: KSP ksp;
150: PetscBool iascii,isstring;
154: if (!viewer) {
155: PetscViewerASCIIGetStdout(((PetscObject)snes)->comm,&viewer);
156: }
160: PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
161: PetscTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
162: if (iascii) {
163: PetscObjectPrintClassNamePrefixType((PetscObject)snes,viewer,"SNES Object");
164: if (snes->ops->view) {
165: PetscViewerASCIIPushTab(viewer);
166: (*snes->ops->view)(snes,viewer);
167: PetscViewerASCIIPopTab(viewer);
168: }
169: PetscViewerASCIIPrintf(viewer," maximum iterations=%D, maximum function evaluations=%D\n",snes->max_its,snes->max_funcs);
170: PetscViewerASCIIPrintf(viewer," tolerances: relative=%G, absolute=%G, solution=%G\n",
171: snes->rtol,snes->abstol,snes->xtol);
172: PetscViewerASCIIPrintf(viewer," total number of linear solver iterations=%D\n",snes->linear_its);
173: PetscViewerASCIIPrintf(viewer," total number of function evaluations=%D\n",snes->nfuncs);
174: if (snes->ksp_ewconv) {
175: kctx = (SNESKSPEW *)snes->kspconvctx;
176: if (kctx) {
177: PetscViewerASCIIPrintf(viewer," Eisenstat-Walker computation of KSP relative tolerance (version %D)\n",kctx->version);
178: PetscViewerASCIIPrintf(viewer," rtol_0=%G, rtol_max=%G, threshold=%G\n",kctx->rtol_0,kctx->rtol_max,kctx->threshold);
179: PetscViewerASCIIPrintf(viewer," gamma=%G, alpha=%G, alpha2=%G\n",kctx->gamma,kctx->alpha,kctx->alpha2);
180: }
181: }
182: if (snes->lagpreconditioner == -1) {
183: PetscViewerASCIIPrintf(viewer," Preconditioned is never rebuilt\n");
184: } else if (snes->lagpreconditioner > 1) {
185: PetscViewerASCIIPrintf(viewer," Preconditioned is rebuilt every %D new Jacobians\n",snes->lagpreconditioner);
186: }
187: if (snes->lagjacobian == -1) {
188: PetscViewerASCIIPrintf(viewer," Jacobian is never rebuilt\n");
189: } else if (snes->lagjacobian > 1) {
190: PetscViewerASCIIPrintf(viewer," Jacobian is rebuilt every %D SNES iterations\n",snes->lagjacobian);
191: }
192: } else if (isstring) {
193: const char *type;
194: SNESGetType(snes,&type);
195: PetscViewerStringSPrintf(viewer," %-3.3s",type);
196: }
197: if (snes->pc && snes->usespc) {
198: PetscViewerASCIIPushTab(viewer);
199: SNESView(snes->pc, viewer);
200: PetscViewerASCIIPopTab(viewer);
201: }
202: if (snes->usesksp) {
203: SNESGetKSP(snes,&ksp);
204: PetscViewerASCIIPushTab(viewer);
205: KSPView(ksp,viewer);
206: PetscViewerASCIIPopTab(viewer);
207: }
208: return(0);
209: }
211: /*
212: We retain a list of functions that also take SNES command
213: line options. These are called at the end SNESSetFromOptions()
214: */
215: #define MAXSETFROMOPTIONS 5
216: static PetscInt numberofsetfromoptions;
217: static PetscErrorCode (*othersetfromoptions[MAXSETFROMOPTIONS])(SNES);
221: /*@C
222: SNESAddOptionsChecker - Adds an additional function to check for SNES options.
224: Not Collective
226: Input Parameter:
227: . snescheck - function that checks for options
229: Level: developer
231: .seealso: SNESSetFromOptions()
232: @*/
233: PetscErrorCode SNESAddOptionsChecker(PetscErrorCode (*snescheck)(SNES))
234: {
236: if (numberofsetfromoptions >= MAXSETFROMOPTIONS) {
237: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Too many options checkers, only %D allowed", MAXSETFROMOPTIONS);
238: }
239: othersetfromoptions[numberofsetfromoptions++] = snescheck;
240: return(0);
241: }
243: extern PetscErrorCode SNESDefaultMatrixFreeCreate2(SNES,Vec,Mat*);
247: static PetscErrorCode SNESSetUpMatrixFree_Private(SNES snes, PetscBool hasOperator, PetscInt version)
248: {
249: Mat J;
250: KSP ksp;
251: PC pc;
252: PetscBool match;
258: if(!snes->vec_func && (snes->jacobian || snes->jacobian_pre)) {
259: Mat A = snes->jacobian, B = snes->jacobian_pre;
260: MatGetVecs(A ? A : B, PETSC_NULL,&snes->vec_func);
261: }
263: if (version == 1) {
264: MatCreateSNESMF(snes,&J);
265: MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);
266: MatSetFromOptions(J);
267: } else if (version == 2) {
268: if (!snes->vec_func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"SNESSetFunction() must be called first");
269: #if !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128)
270: SNESDefaultMatrixFreeCreate2(snes,snes->vec_func,&J);
271: #else
272: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "matrix-free operator rutines (version 2)");
273: #endif
274: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "matrix-free operator rutines, only version 1 and 2");
275:
276: PetscInfo1(snes,"Setting default matrix-free operator routines (version %D)\n", version);
277: if (hasOperator) {
278: /* This version replaces the user provided Jacobian matrix with a
279: matrix-free version but still employs the user-provided preconditioner matrix. */
280: SNESSetJacobian(snes,J,0,0,0);
281: } else {
282: /* This version replaces both the user-provided Jacobian and the user-
283: provided preconditioner matrix with the default matrix free version. */
284: SNESSetJacobian(snes,J,J,MatMFFDComputeJacobian,snes->funP);
285: /* Force no preconditioner */
286: SNESGetKSP(snes,&ksp);
287: KSPGetPC(ksp,&pc);
288: PetscTypeCompare((PetscObject)pc,PCSHELL,&match);
289: if (!match) {
290: PetscInfo(snes,"Setting default matrix-free preconditioner routines\nThat is no preconditioner is being used\n");
291: PCSetType(pc,PCNONE);
292: }
293: }
294: MatDestroy(&J);
295: return(0);
296: }
300: /*@
301: SNESSetFromOptions - Sets various SNES and KSP parameters from user options.
303: Collective on SNES
305: Input Parameter:
306: . snes - the SNES context
308: Options Database Keys:
309: + -snes_type <type> - ls, tr, ngmres, ncg, richardson, qn, vi, fas
310: . -snes_stol - convergence tolerance in terms of the norm
311: of the change in the solution between steps
312: . -snes_atol <abstol> - absolute tolerance of residual norm
313: . -snes_rtol <rtol> - relative decrease in tolerance norm from initial
314: . -snes_max_it <max_it> - maximum number of iterations
315: . -snes_max_funcs <max_funcs> - maximum number of function evaluations
316: . -snes_max_fail <max_fail> - maximum number of line search failures allowed before stopping, default is none
317: . -snes_max_linear_solve_fail - number of linear solver failures before SNESSolve() stops
318: . -snes_lag_preconditioner <lag> - how often preconditioner is rebuilt (use -1 to never rebuild)
319: . -snes_lag_jacobian <lag> - how often Jacobian is rebuilt (use -1 to never rebuild)
320: . -snes_trtol <trtol> - trust region tolerance
321: . -snes_no_convergence_test - skip convergence test in nonlinear
322: solver; hence iterations will continue until max_it
323: or some other criterion is reached. Saves expense
324: of convergence test
325: . -snes_monitor <optional filename> - prints residual norm at each iteration. if no
326: filename given prints to stdout
327: . -snes_monitor_solution - plots solution at each iteration
328: . -snes_monitor_residual - plots residual (not its norm) at each iteration
329: . -snes_monitor_solution_update - plots update to solution at each iteration
330: . -snes_monitor_draw - plots residual norm at each iteration
331: . -snes_fd - use finite differences to compute Jacobian; very slow, only for testing
332: . -snes_mf_ksp_monitor - if using matrix-free multiply then print h at each KSP iteration
333: - -snes_converged_reason - print the reason for convergence/divergence after each solve
335: Options Database for Eisenstat-Walker method:
336: + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence
337: . -snes_ksp_ew_version ver - version of Eisenstat-Walker method
338: . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0
339: . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax
340: . -snes_ksp_ew_gamma <gamma> - Sets gamma
341: . -snes_ksp_ew_alpha <alpha> - Sets alpha
342: . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2
343: - -snes_ksp_ew_threshold <threshold> - Sets threshold
345: Notes:
346: To see all options, run your program with the -help option or consult
347: the <A href="../../docs/manual.pdf#nameddest=ch_snes">SNES chapter of the users manual</A>.
349: Level: beginner
351: .keywords: SNES, nonlinear, set, options, database
353: .seealso: SNESSetOptionsPrefix()
354: @*/
355: PetscErrorCode SNESSetFromOptions(SNES snes)
356: {
357: PetscBool flg,set,mf,mf_operator,pcset;
358: PetscInt i,indx,lag,mf_version,grids;
359: MatStructure matflag;
360: const char *deft = SNESLS;
361: const char *convtests[] = {"default","skip"};
362: SNESKSPEW *kctx = NULL;
363: char type[256], monfilename[PETSC_MAX_PATH_LEN];
364: const char *optionsprefix;
365: PetscViewer monviewer;
366: PetscErrorCode ierr;
371: if (!SNESRegisterAllCalled) {SNESRegisterAll(PETSC_NULL);}
372: PetscObjectOptionsBegin((PetscObject)snes);
373: if (((PetscObject)snes)->type_name) { deft = ((PetscObject)snes)->type_name; }
374: PetscOptionsList("-snes_type","Nonlinear solver method","SNESSetType",SNESList,deft,type,256,&flg);
375: if (flg) {
376: SNESSetType(snes,type);
377: } else if (!((PetscObject)snes)->type_name) {
378: SNESSetType(snes,deft);
379: }
380: /* not used here, but called so will go into help messaage */
381: PetscOptionsName("-snes_view","Print detailed information on solver used","SNESView",0);
383: PetscOptionsReal("-snes_stol","Stop if step length less than","SNESSetTolerances",snes->xtol,&snes->xtol,0);
384: PetscOptionsReal("-snes_atol","Stop if function norm less than","SNESSetTolerances",snes->abstol,&snes->abstol,0);
386: PetscOptionsReal("-snes_rtol","Stop if decrease in function norm less than","SNESSetTolerances",snes->rtol,&snes->rtol,0);
387: PetscOptionsInt("-snes_max_it","Maximum iterations","SNESSetTolerances",snes->max_its,&snes->max_its,PETSC_NULL);
388: PetscOptionsInt("-snes_max_funcs","Maximum function evaluations","SNESSetTolerances",snes->max_funcs,&snes->max_funcs,PETSC_NULL);
389: PetscOptionsInt("-snes_max_fail","Maximum failures","SNESSetTolerances",snes->maxFailures,&snes->maxFailures,PETSC_NULL);
390: PetscOptionsInt("-snes_max_linear_solve_fail","Maximum failures in linear solves allowed","SNESSetMaxLinearSolveFailures",snes->maxLinearSolveFailures,&snes->maxLinearSolveFailures,PETSC_NULL);
391: PetscOptionsBool("-snes_error_if_not_converged","Generate error if solver does not converge","SNESSetErrorIfNotConverged",snes->errorifnotconverged,&snes->errorifnotconverged,PETSC_NULL);
393: PetscOptionsInt("-snes_lag_preconditioner","How often to rebuild preconditioner","SNESSetLagPreconditioner",snes->lagpreconditioner,&lag,&flg);
394: if (flg) {
395: SNESSetLagPreconditioner(snes,lag);
396: }
397: PetscOptionsInt("-snes_lag_jacobian","How often to rebuild Jacobian","SNESSetLagJacobian",snes->lagjacobian,&lag,&flg);
398: if (flg) {
399: SNESSetLagJacobian(snes,lag);
400: }
401: PetscOptionsInt("-snes_grid_sequence","Use grid sequencing to generate initial guess","SNESSetGridSequence",snes->gridsequence,&grids,&flg);
402: if (flg) {
403: SNESSetGridSequence(snes,grids);
404: }
406: PetscOptionsEList("-snes_convergence_test","Convergence test","SNESSetConvergenceTest",convtests,2,"default",&indx,&flg);
407: if (flg) {
408: switch (indx) {
409: case 0: SNESSetConvergenceTest(snes,SNESDefaultConverged,PETSC_NULL,PETSC_NULL); break;
410: case 1: SNESSetConvergenceTest(snes,SNESSkipConverged,PETSC_NULL,PETSC_NULL); break;
411: }
412: }
414: PetscOptionsBool("-snes_converged_reason","Print reason for converged or diverged","SNESSolve",snes->printreason,&snes->printreason,PETSC_NULL);
416: kctx = (SNESKSPEW *)snes->kspconvctx;
418: PetscOptionsBool("-snes_ksp_ew","Use Eisentat-Walker linear system convergence test","SNESKSPSetUseEW",snes->ksp_ewconv,&snes->ksp_ewconv,PETSC_NULL);
420: PetscOptionsInt("-snes_ksp_ew_version","Version 1, 2 or 3","SNESKSPSetParametersEW",kctx->version,&kctx->version,0);
421: PetscOptionsReal("-snes_ksp_ew_rtol0","0 <= rtol0 < 1","SNESKSPSetParametersEW",kctx->rtol_0,&kctx->rtol_0,0);
422: PetscOptionsReal("-snes_ksp_ew_rtolmax","0 <= rtolmax < 1","SNESKSPSetParametersEW",kctx->rtol_max,&kctx->rtol_max,0);
423: PetscOptionsReal("-snes_ksp_ew_gamma","0 <= gamma <= 1","SNESKSPSetParametersEW",kctx->gamma,&kctx->gamma,0);
424: PetscOptionsReal("-snes_ksp_ew_alpha","1 < alpha <= 2","SNESKSPSetParametersEW",kctx->alpha,&kctx->alpha,0);
425: PetscOptionsReal("-snes_ksp_ew_alpha2","alpha2","SNESKSPSetParametersEW",kctx->alpha2,&kctx->alpha2,0);
426: PetscOptionsReal("-snes_ksp_ew_threshold","0 < threshold < 1","SNESKSPSetParametersEW",kctx->threshold,&kctx->threshold,0);
428: flg = PETSC_FALSE;
429: PetscOptionsBool("-snes_monitor_cancel","Remove all monitors","SNESMonitorCancel",flg,&flg,PETSC_NULL);
430: if (flg) {SNESMonitorCancel(snes);}
432: PetscOptionsString("-snes_monitor","Monitor norm of function","SNESMonitorSet","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);
433: if (flg) {
434: PetscViewerASCIIOpen(((PetscObject)snes)->comm,monfilename,&monviewer);
435: SNESMonitorSet(snes,SNESMonitorDefault,monviewer,(PetscErrorCode (*)(void**))PetscViewerDestroy);
436: }
438: PetscOptionsString("-snes_monitor_range","Monitor range of elements of function","SNESMonitorSet","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);
439: if (flg) {
440: SNESMonitorSet(snes,SNESMonitorRange,0,0);
441: }
443: PetscOptionsString("-snes_ratiomonitor","Monitor ratios of norms of function","SNESMonitorSetRatio","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);
444: if (flg) {
445: PetscViewerASCIIOpen(((PetscObject)snes)->comm,monfilename,&monviewer);
446: SNESMonitorSetRatio(snes,monviewer);
447: }
449: PetscOptionsString("-snes_monitor_short","Monitor norm of function (fewer digits)","SNESMonitorSet","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);
450: if (flg) {
451: PetscViewerASCIIOpen(((PetscObject)snes)->comm,monfilename,&monviewer);
452: SNESMonitorSet(snes,SNESMonitorDefaultShort,monviewer,(PetscErrorCode (*)(void**))PetscViewerDestroy);
453: }
455: PetscOptionsString("-snes_monitor_python","Use Python function","SNESMonitorSet",0,monfilename,PETSC_MAX_PATH_LEN,&flg);
456: if (flg) {PetscPythonMonitorSet((PetscObject)snes,monfilename);}
458: flg = PETSC_FALSE;
459: PetscOptionsBool("-snes_monitor_solution","Plot solution at each iteration","SNESMonitorSolution",flg,&flg,PETSC_NULL);
460: if (flg) {SNESMonitorSet(snes,SNESMonitorSolution,0,0);}
461: flg = PETSC_FALSE;
462: PetscOptionsBool("-snes_monitor_solution_update","Plot correction at each iteration","SNESMonitorSolutionUpdate",flg,&flg,PETSC_NULL);
463: if (flg) {SNESMonitorSet(snes,SNESMonitorSolutionUpdate,0,0);}
464: flg = PETSC_FALSE;
465: PetscOptionsBool("-snes_monitor_residual","Plot residual at each iteration","SNESMonitorResidual",flg,&flg,PETSC_NULL);
466: if (flg) {SNESMonitorSet(snes,SNESMonitorResidual,0,0);}
467: flg = PETSC_FALSE;
468: PetscOptionsBool("-snes_monitor_draw","Plot function norm at each iteration","SNESMonitorLG",flg,&flg,PETSC_NULL);
469: if (flg) {SNESMonitorSet(snes,SNESMonitorLG,PETSC_NULL,PETSC_NULL);}
470: flg = PETSC_FALSE;
471: PetscOptionsBool("-snes_monitor_range_draw","Plot function range at each iteration","SNESMonitorLG",flg,&flg,PETSC_NULL);
472: if (flg) {SNESMonitorSet(snes,SNESMonitorLGRange,PETSC_NULL,PETSC_NULL);}
474: flg = PETSC_FALSE;
475: PetscOptionsBool("-snes_fd","Use finite differences (slow) to compute Jacobian","SNESDefaultComputeJacobian",flg,&flg,PETSC_NULL);
476: if (flg) {
477: SNESSetJacobian(snes,snes->jacobian,snes->jacobian_pre,SNESDefaultComputeJacobian,snes->funP);
478: PetscInfo(snes,"Setting default finite difference Jacobian matrix\n");
479: }
481: mf = mf_operator = PETSC_FALSE;
482: flg = PETSC_FALSE;
483: PetscOptionsBool("-snes_mf_operator","Use a Matrix-Free Jacobian with user-provided preconditioner matrix","MatCreateSNESMF",PETSC_FALSE,&mf_operator,&flg);
484: if (flg && mf_operator) {
485: snes->mf_operator = PETSC_TRUE;
486: mf = PETSC_TRUE;
487: }
488: flg = PETSC_FALSE;
489: PetscOptionsBool("-snes_mf","Use a Matrix-Free Jacobian with no preconditioner matrix","MatCreateSNESMF",PETSC_FALSE,&mf,&flg);
490: if (!flg && mf_operator) mf = PETSC_TRUE;
491: mf_version = 1;
492: PetscOptionsInt("-snes_mf_version","Matrix-Free routines version 1 or 2","None",mf_version,&mf_version,0);
495: /* GS Options */
496: PetscOptionsInt("-snes_gs_sweeps","Number of sweeps of GS to apply","SNESComputeGS",snes->gssweeps,&snes->gssweeps,PETSC_NULL);
498: /* line search options */
499: PetscOptionsReal("-snes_ls_alpha","Constant function norm must decrease by","None",snes->ls_alpha,&snes->ls_alpha,0);
500: PetscOptionsReal("-snes_ls_maxstep","Step must be less than","None",snes->maxstep,&snes->maxstep,0);
501: PetscOptionsReal("-snes_ls_steptol","Minimum lambda allowed","None",snes->steptol,&snes->steptol,0);
502: PetscOptionsReal("-snes_ls_damping","Damping parameter","SNES",snes->damping,&snes->damping,&flg);
503: PetscOptionsInt("-snes_ls_it" ,"Line search iterations","SNES",snes->ls_its,&snes->ls_its,&flg);
504: PetscOptionsBool("-snes_ls_monitor","Print progress of line searches","SNESLineSearchSetMonitor",snes->ls_monitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);
505: if (set) {SNESLineSearchSetMonitor(snes,flg);}
506: PetscOptionsEnum("-snes_ls","Line search used","SNESLineSearchSet",SNESLineSearchTypes,(PetscEnum)snes->ls_type,(PetscEnum*)&indx,&flg);
507: if (flg) {
508: SNESLineSearchSetType(snes,(SNESLineSearchType)indx);
509: }
510: flg = snes->ops->precheckstep == SNESLineSearchPreCheckPicard ? PETSC_TRUE : PETSC_FALSE;
511: PetscOptionsBool("-snes_ls_precheck_picard","Use a correction that sometimes improves convergence of Picard iteration","SNESLineSearchPreCheckPicard",flg,&flg,&set);
512: if (set) {
513: if (flg) {
514: snes->precheck_picard_angle = 10.; /* correction only active if angle is less than 10 degrees */
515: PetscOptionsReal("-snes_ls_precheck_picard_angle","Maximum angle at which to activate the correction","none",snes->precheck_picard_angle,&snes->precheck_picard_angle,PETSC_NULL);
516: SNESLineSearchSetPreCheck(snes,SNESLineSearchPreCheckPicard,&snes->precheck_picard_angle);
517: } else {
518: SNESLineSearchSetPreCheck(snes,PETSC_NULL,PETSC_NULL);
519: }
520: }
522: for(i = 0; i < numberofsetfromoptions; i++) {
523: (*othersetfromoptions[i])(snes);
524: }
526: if (snes->ops->setfromoptions) {
527: (*snes->ops->setfromoptions)(snes);
528: }
530: /* process any options handlers added with PetscObjectAddOptionsHandler() */
531: PetscObjectProcessOptionsHandlers((PetscObject)snes);
532: PetscOptionsEnd();
534: if (mf) { SNESSetUpMatrixFree_Private(snes, mf_operator, mf_version); }
536: if (!snes->ksp) {SNESGetKSP(snes,&snes->ksp);}
537: KSPGetOperators(snes->ksp,PETSC_NULL,PETSC_NULL,&matflag);
538: KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,matflag);
539: KSPSetFromOptions(snes->ksp);
541: /* if someone has set the SNES PC type, create it. */
542: SNESGetOptionsPrefix(snes, &optionsprefix);
543: PetscOptionsHasName(optionsprefix, "-npc_snes_type", &pcset);
544: if (pcset && (!snes->pc)) {
545: SNESGetPC(snes, &snes->pc);
546: }
547: if (snes->pc) {
548: SNESSetOptionsPrefix(snes->pc, optionsprefix);
549: SNESAppendOptionsPrefix(snes->pc, "npc_");
550: SNESSetDM(snes->pc, snes->dm);
551: SNESSetGS(snes->pc, snes->ops->computegs, snes->gsP);
552: /* Should we make a duplicate vector and matrix? Leave the DM to make it? */
553: SNESSetFunction(snes->pc, snes->vec_func, snes->ops->computefunction, snes->funP);
554: SNESSetJacobian(snes->pc, snes->jacobian, snes->jacobian_pre, snes->ops->computejacobian, snes->jacP);
555: SNESSetFromOptions(snes->pc);
556: }
557: return(0);
558: }
562: /*@
563: SNESSetComputeApplicationContext - Sets an optional function to compute a user-defined context for
564: the nonlinear solvers.
566: Logically Collective on SNES
568: Input Parameters:
569: + snes - the SNES context
570: . compute - function to compute the context
571: - destroy - function to destroy the context
573: Level: intermediate
575: .keywords: SNES, nonlinear, set, application, context
577: .seealso: SNESGetApplicationContext(), SNESSetComputeApplicationContext(), SNESGetApplicationContext()
578: @*/
579: PetscErrorCode SNESSetComputeApplicationContext(SNES snes,PetscErrorCode (*compute)(SNES,void**),PetscErrorCode (*destroy)(void**))
580: {
583: snes->ops->usercompute = compute;
584: snes->ops->userdestroy = destroy;
585: return(0);
586: }
590: /*@
591: SNESSetApplicationContext - Sets the optional user-defined context for
592: the nonlinear solvers.
594: Logically Collective on SNES
596: Input Parameters:
597: + snes - the SNES context
598: - usrP - optional user context
600: Level: intermediate
602: .keywords: SNES, nonlinear, set, application, context
604: .seealso: SNESGetApplicationContext(), SNESSetApplicationContext()
605: @*/
606: PetscErrorCode SNESSetApplicationContext(SNES snes,void *usrP)
607: {
609: KSP ksp;
613: SNESGetKSP(snes,&ksp);
614: KSPSetApplicationContext(ksp,usrP);
615: snes->user = usrP;
616: return(0);
617: }
621: /*@
622: SNESGetApplicationContext - Gets the user-defined context for the
623: nonlinear solvers.
625: Not Collective
627: Input Parameter:
628: . snes - SNES context
630: Output Parameter:
631: . usrP - user context
633: Level: intermediate
635: .keywords: SNES, nonlinear, get, application, context
637: .seealso: SNESSetApplicationContext()
638: @*/
639: PetscErrorCode SNESGetApplicationContext(SNES snes,void *usrP)
640: {
643: *(void**)usrP = snes->user;
644: return(0);
645: }
649: /*@
650: SNESGetIterationNumber - Gets the number of nonlinear iterations completed
651: at this time.
653: Not Collective
655: Input Parameter:
656: . snes - SNES context
658: Output Parameter:
659: . iter - iteration number
661: Notes:
662: For example, during the computation of iteration 2 this would return 1.
664: This is useful for using lagged Jacobians (where one does not recompute the
665: Jacobian at each SNES iteration). For example, the code
666: .vb
667: SNESGetIterationNumber(snes,&it);
668: if (!(it % 2)) {
669: [compute Jacobian here]
670: }
671: .ve
672: can be used in your ComputeJacobian() function to cause the Jacobian to be
673: recomputed every second SNES iteration.
675: Level: intermediate
677: .keywords: SNES, nonlinear, get, iteration, number,
679: .seealso: SNESGetFunctionNorm(), SNESGetLinearSolveIterations()
680: @*/
681: PetscErrorCode SNESGetIterationNumber(SNES snes,PetscInt* iter)
682: {
686: *iter = snes->iter;
687: return(0);
688: }
692: /*@
693: SNESSetIterationNumber - Sets the current iteration number.
695: Not Collective
697: Input Parameter:
698: . snes - SNES context
699: . iter - iteration number
701: Level: developer
703: .keywords: SNES, nonlinear, set, iteration, number,
705: .seealso: SNESGetFunctionNorm(), SNESGetLinearSolveIterations()
706: @*/
707: PetscErrorCode SNESSetIterationNumber(SNES snes,PetscInt iter)
708: {
713: PetscObjectTakeAccess(snes);
714: snes->iter = iter;
715: PetscObjectGrantAccess(snes);
716: return(0);
717: }
721: /*@
722: SNESGetFunctionNorm - Gets the norm of the current function that was set
723: with SNESSSetFunction().
725: Collective on SNES
727: Input Parameter:
728: . snes - SNES context
730: Output Parameter:
731: . fnorm - 2-norm of function
733: Level: intermediate
735: .keywords: SNES, nonlinear, get, function, norm
737: .seealso: SNESGetFunction(), SNESGetIterationNumber(), SNESGetLinearSolveIterations()
738: @*/
739: PetscErrorCode SNESGetFunctionNorm(SNES snes,PetscReal *fnorm)
740: {
744: *fnorm = snes->norm;
745: return(0);
746: }
751: /*@
752: SNESSetFunctionNorm - Sets the 2-norm of the current function computed using VecNorm().
754: Collective on SNES
756: Input Parameter:
757: . snes - SNES context
758: . fnorm - 2-norm of function
760: Level: developer
762: .keywords: SNES, nonlinear, set, function, norm
764: .seealso: SNESSetFunction(), SNESSetIterationNumber(), VecNorm().
765: @*/
766: PetscErrorCode SNESSetFunctionNorm(SNES snes,PetscReal fnorm)
767: {
773: PetscObjectTakeAccess(snes);
774: snes->norm = fnorm;
775: PetscObjectGrantAccess(snes);
776: return(0);
777: }
781: /*@
782: SNESGetNonlinearStepFailures - Gets the number of unsuccessful steps
783: attempted by the nonlinear solver.
785: Not Collective
787: Input Parameter:
788: . snes - SNES context
790: Output Parameter:
791: . nfails - number of unsuccessful steps attempted
793: Notes:
794: This counter is reset to zero for each successive call to SNESSolve().
796: Level: intermediate
798: .keywords: SNES, nonlinear, get, number, unsuccessful, steps
800: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures(),
801: SNESSetMaxNonlinearStepFailures(), SNESGetMaxNonlinearStepFailures()
802: @*/
803: PetscErrorCode SNESGetNonlinearStepFailures(SNES snes,PetscInt* nfails)
804: {
808: *nfails = snes->numFailures;
809: return(0);
810: }
814: /*@
815: SNESSetMaxNonlinearStepFailures - Sets the maximum number of unsuccessful steps
816: attempted by the nonlinear solver before it gives up.
818: Not Collective
820: Input Parameters:
821: + snes - SNES context
822: - maxFails - maximum of unsuccessful steps
824: Level: intermediate
826: .keywords: SNES, nonlinear, set, maximum, unsuccessful, steps
828: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures(),
829: SNESGetMaxNonlinearStepFailures(), SNESGetNonlinearStepFailures()
830: @*/
831: PetscErrorCode SNESSetMaxNonlinearStepFailures(SNES snes, PetscInt maxFails)
832: {
835: snes->maxFailures = maxFails;
836: return(0);
837: }
841: /*@
842: SNESGetMaxNonlinearStepFailures - Gets the maximum number of unsuccessful steps
843: attempted by the nonlinear solver before it gives up.
845: Not Collective
847: Input Parameter:
848: . snes - SNES context
850: Output Parameter:
851: . maxFails - maximum of unsuccessful steps
853: Level: intermediate
855: .keywords: SNES, nonlinear, get, maximum, unsuccessful, steps
857: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures(),
858: SNESSetMaxNonlinearStepFailures(), SNESGetNonlinearStepFailures()
859:
860: @*/
861: PetscErrorCode SNESGetMaxNonlinearStepFailures(SNES snes, PetscInt *maxFails)
862: {
866: *maxFails = snes->maxFailures;
867: return(0);
868: }
872: /*@
873: SNESGetNumberFunctionEvals - Gets the number of user provided function evaluations
874: done by SNES.
876: Not Collective
878: Input Parameter:
879: . snes - SNES context
881: Output Parameter:
882: . nfuncs - number of evaluations
884: Level: intermediate
886: .keywords: SNES, nonlinear, get, maximum, unsuccessful, steps
888: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures()
889: @*/
890: PetscErrorCode SNESGetNumberFunctionEvals(SNES snes, PetscInt *nfuncs)
891: {
895: *nfuncs = snes->nfuncs;
896: return(0);
897: }
901: /*@
902: SNESGetLinearSolveFailures - Gets the number of failed (non-converged)
903: linear solvers.
905: Not Collective
907: Input Parameter:
908: . snes - SNES context
910: Output Parameter:
911: . nfails - number of failed solves
913: Notes:
914: This counter is reset to zero for each successive call to SNESSolve().
916: Level: intermediate
918: .keywords: SNES, nonlinear, get, number, unsuccessful, steps
920: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures()
921: @*/
922: PetscErrorCode SNESGetLinearSolveFailures(SNES snes,PetscInt* nfails)
923: {
927: *nfails = snes->numLinearSolveFailures;
928: return(0);
929: }
933: /*@
934: SNESSetMaxLinearSolveFailures - the number of failed linear solve attempts
935: allowed before SNES returns with a diverged reason of SNES_DIVERGED_LINEAR_SOLVE
937: Logically Collective on SNES
939: Input Parameters:
940: + snes - SNES context
941: - maxFails - maximum allowed linear solve failures
943: Level: intermediate
945: Notes: By default this is 0; that is SNES returns on the first failed linear solve
947: .keywords: SNES, nonlinear, set, maximum, unsuccessful, steps
949: .seealso: SNESGetLinearSolveFailures(), SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations()
950: @*/
951: PetscErrorCode SNESSetMaxLinearSolveFailures(SNES snes, PetscInt maxFails)
952: {
956: snes->maxLinearSolveFailures = maxFails;
957: return(0);
958: }
962: /*@
963: SNESGetMaxLinearSolveFailures - gets the maximum number of linear solve failures that
964: are allowed before SNES terminates
966: Not Collective
968: Input Parameter:
969: . snes - SNES context
971: Output Parameter:
972: . maxFails - maximum of unsuccessful solves allowed
974: Level: intermediate
976: Notes: By default this is 1; that is SNES returns on the first failed linear solve
978: .keywords: SNES, nonlinear, get, maximum, unsuccessful, steps
980: .seealso: SNESGetLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(),
981: @*/
982: PetscErrorCode SNESGetMaxLinearSolveFailures(SNES snes, PetscInt *maxFails)
983: {
987: *maxFails = snes->maxLinearSolveFailures;
988: return(0);
989: }
993: /*@
994: SNESGetLinearSolveIterations - Gets the total number of linear iterations
995: used by the nonlinear solver.
997: Not Collective
999: Input Parameter:
1000: . snes - SNES context
1002: Output Parameter:
1003: . lits - number of linear iterations
1005: Notes:
1006: This counter is reset to zero for each successive call to SNESSolve().
1008: Level: intermediate
1010: .keywords: SNES, nonlinear, get, number, linear, iterations
1012: .seealso: SNESGetIterationNumber(), SNESGetFunctionNorm(), SNESGetLinearSolveFailures(), SNESGetMaxLinearSolveFailures()
1013: @*/
1014: PetscErrorCode SNESGetLinearSolveIterations(SNES snes,PetscInt* lits)
1015: {
1019: *lits = snes->linear_its;
1020: return(0);
1021: }
1025: /*@
1026: SNESGetKSP - Returns the KSP context for a SNES solver.
1028: Not Collective, but if SNES object is parallel, then KSP object is parallel
1030: Input Parameter:
1031: . snes - the SNES context
1033: Output Parameter:
1034: . ksp - the KSP context
1036: Notes:
1037: The user can then directly manipulate the KSP context to set various
1038: options, etc. Likewise, the user can then extract and manipulate the
1039: PC contexts as well.
1041: Level: beginner
1043: .keywords: SNES, nonlinear, get, KSP, context
1045: .seealso: KSPGetPC(), SNESCreate(), KSPCreate(), SNESSetKSP()
1046: @*/
1047: PetscErrorCode SNESGetKSP(SNES snes,KSP *ksp)
1048: {
1055: if (!snes->ksp) {
1056: KSPCreate(((PetscObject)snes)->comm,&snes->ksp);
1057: PetscObjectIncrementTabLevel((PetscObject)snes->ksp,(PetscObject)snes,1);
1058: PetscLogObjectParent(snes,snes->ksp);
1059: }
1060: *ksp = snes->ksp;
1061: return(0);
1062: }
1066: /*@
1067: SNESSetKSP - Sets a KSP context for the SNES object to use
1069: Not Collective, but the SNES and KSP objects must live on the same MPI_Comm
1071: Input Parameters:
1072: + snes - the SNES context
1073: - ksp - the KSP context
1075: Notes:
1076: The SNES object already has its KSP object, you can obtain with SNESGetKSP()
1077: so this routine is rarely needed.
1079: The KSP object that is already in the SNES object has its reference count
1080: decreased by one.
1082: Level: developer
1084: .keywords: SNES, nonlinear, get, KSP, context
1086: .seealso: KSPGetPC(), SNESCreate(), KSPCreate(), SNESSetKSP()
1087: @*/
1088: PetscErrorCode SNESSetKSP(SNES snes,KSP ksp)
1089: {
1096: PetscObjectReference((PetscObject)ksp);
1097: if (snes->ksp) {PetscObjectDereference((PetscObject)snes->ksp);}
1098: snes->ksp = ksp;
1099: return(0);
1100: }
1102: #if 0
1105: static PetscErrorCode SNESPublish_Petsc(PetscObject obj)
1106: {
1108: return(0);
1109: }
1110: #endif
1112: /* -----------------------------------------------------------*/
1115: /*@
1116: SNESCreate - Creates a nonlinear solver context.
1118: Collective on MPI_Comm
1120: Input Parameters:
1121: . comm - MPI communicator
1123: Output Parameter:
1124: . outsnes - the new SNES context
1126: Options Database Keys:
1127: + -snes_mf - Activates default matrix-free Jacobian-vector products,
1128: and no preconditioning matrix
1129: . -snes_mf_operator - Activates default matrix-free Jacobian-vector
1130: products, and a user-provided preconditioning matrix
1131: as set by SNESSetJacobian()
1132: - -snes_fd - Uses (slow!) finite differences to compute Jacobian
1134: Level: beginner
1136: .keywords: SNES, nonlinear, create, context
1138: .seealso: SNESSolve(), SNESDestroy(), SNES, SNESSetLagPreconditioner()
1140: @*/
1141: PetscErrorCode SNESCreate(MPI_Comm comm,SNES *outsnes)
1142: {
1143: PetscErrorCode ierr;
1144: SNES snes;
1145: SNESKSPEW *kctx;
1149: *outsnes = PETSC_NULL;
1150: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
1151: SNESInitializePackage(PETSC_NULL);
1152: #endif
1154: PetscHeaderCreate(snes,_p_SNES,struct _SNESOps,SNES_CLASSID,0,"SNES","Nonlinear solver","SNES",comm,SNESDestroy,SNESView);
1156: snes->ops->converged = SNESDefaultConverged;
1157: snes->usesksp = PETSC_TRUE;
1158: snes->max_its = 50;
1159: snes->max_funcs = 10000;
1160: snes->norm = 0.0;
1161: snes->rtol = 1.e-8;
1162: snes->ttol = 0.0;
1163: snes->abstol = 1.e-50;
1164: snes->xtol = 1.e-8;
1165: snes->deltatol = 1.e-12;
1166: snes->nfuncs = 0;
1167: snes->numFailures = 0;
1168: snes->maxFailures = 1;
1169: snes->linear_its = 0;
1170: snes->lagjacobian = 1;
1171: snes->lagpreconditioner = 1;
1172: snes->numbermonitors = 0;
1173: snes->data = 0;
1174: snes->setupcalled = PETSC_FALSE;
1175: snes->ksp_ewconv = PETSC_FALSE;
1176: snes->nwork = 0;
1177: snes->work = 0;
1178: snes->nvwork = 0;
1179: snes->vwork = 0;
1180: snes->conv_hist_len = 0;
1181: snes->conv_hist_max = 0;
1182: snes->conv_hist = PETSC_NULL;
1183: snes->conv_hist_its = PETSC_NULL;
1184: snes->conv_hist_reset = PETSC_TRUE;
1185: snes->reason = SNES_CONVERGED_ITERATING;
1186: snes->gssweeps = 1;
1188: /* initialize the line search options */
1189: snes->ls_type = SNES_LS_BASIC;
1190: snes->ls_its = 1;
1191: snes->damping = 1.0;
1192: snes->maxstep = 1e8;
1193: snes->steptol = 1e-12;
1194: snes->ls_alpha = 1e-4;
1195: snes->ls_monitor = PETSC_NULL;
1197: snes->ops->linesearch = PETSC_NULL;
1198: snes->precheck = PETSC_NULL;
1199: snes->ops->precheckstep = PETSC_NULL;
1200: snes->postcheck = PETSC_NULL;
1201: snes->ops->postcheckstep= PETSC_NULL;
1203: snes->numLinearSolveFailures = 0;
1204: snes->maxLinearSolveFailures = 1;
1206: /* Create context to compute Eisenstat-Walker relative tolerance for KSP */
1207: PetscNewLog(snes,SNESKSPEW,&kctx);
1208: snes->kspconvctx = (void*)kctx;
1209: kctx->version = 2;
1210: kctx->rtol_0 = .3; /* Eisenstat and Walker suggest rtol_0=.5, but
1211: this was too large for some test cases */
1212: kctx->rtol_last = 0.0;
1213: kctx->rtol_max = .9;
1214: kctx->gamma = 1.0;
1215: kctx->alpha = .5*(1.0 + PetscSqrtReal(5.0));
1216: kctx->alpha2 = kctx->alpha;
1217: kctx->threshold = .1;
1218: kctx->lresid_last = 0.0;
1219: kctx->norm_last = 0.0;
1221: *outsnes = snes;
1222: return(0);
1223: }
1227: /*@C
1228: SNESSetFunction - Sets the function evaluation routine and function
1229: vector for use by the SNES routines in solving systems of nonlinear
1230: equations.
1232: Logically Collective on SNES
1234: Input Parameters:
1235: + snes - the SNES context
1236: . r - vector to store function value
1237: . func - function evaluation routine
1238: - ctx - [optional] user-defined context for private data for the
1239: function evaluation routine (may be PETSC_NULL)
1241: Calling sequence of func:
1242: $ func (SNES snes,Vec x,Vec f,void *ctx);
1244: . f - function vector
1245: - ctx - optional user-defined function context
1247: Notes:
1248: The Newton-like methods typically solve linear systems of the form
1249: $ f'(x) x = -f(x),
1250: where f'(x) denotes the Jacobian matrix and f(x) is the function.
1252: Level: beginner
1254: .keywords: SNES, nonlinear, set, function
1256: .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetPicard()
1257: @*/
1258: PetscErrorCode SNESSetFunction(SNES snes,Vec r,PetscErrorCode (*func)(SNES,Vec,Vec,void*),void *ctx)
1259: {
1263: if (r) {
1266: PetscObjectReference((PetscObject)r);
1267: VecDestroy(&snes->vec_func);
1268: snes->vec_func = r;
1269: } else if (!snes->vec_func && snes->dm) {
1270: DMCreateGlobalVector(snes->dm,&snes->vec_func);
1271: }
1272: if (func) snes->ops->computefunction = func;
1273: if (ctx) snes->funP = ctx;
1274: return(0);
1275: }
1280: /*@C
1281: SNESSetGS - Sets the user nonlinear Gauss-Seidel routine for
1282: use with composed nonlinear solvers.
1284: Input Parameters:
1285: + snes - the SNES context
1286: . gsfunc - function evaluation routine
1287: - ctx - [optional] user-defined context for private data for the
1288: smoother evaluation routine (may be PETSC_NULL)
1290: Calling sequence of func:
1291: $ func (SNES snes,Vec x,Vec b,void *ctx);
1293: + X - solution vector
1294: . B - RHS vector
1295: - ctx - optional user-defined Gauss-Seidel context
1297: Notes:
1298: The GS routines are used by the composed nonlinear solver to generate
1299: a problem appropriate update to the solution, particularly FAS.
1301: Level: intermediate
1303: .keywords: SNES, nonlinear, set, Gauss-Seidel
1305: .seealso: SNESGetFunction(), SNESComputeGS()
1306: @*/
1307: PetscErrorCode SNESSetGS(SNES snes, PetscErrorCode (*gsfunc)(SNES,Vec,Vec,void *), void * ctx) {
1309: if (gsfunc) snes->ops->computegs = gsfunc;
1310: if (ctx) snes->gsP = ctx;
1311: return(0);
1312: }
1316: /*@
1317: SNESSetGSSweeps - Sets the number of sweeps of GS to use.
1319: Input Parameters:
1320: + snes - the SNES context
1321: - sweeps - the number of sweeps of GS to perform.
1323: Level: intermediate
1325: .keywords: SNES, nonlinear, set, Gauss-Siedel
1327: .seealso: SNESSetGS(), SNESGetGS(), SNESSetPC(), SNESGetGSSweeps()
1328: @*/
1330: PetscErrorCode SNESSetGSSweeps(SNES snes, PetscInt sweeps) {
1332: snes->gssweeps = sweeps;
1333: return(0);
1334: }
1339: /*@
1340: SNESGetGSSweeps - Gets the number of sweeps GS will use.
1342: Input Parameters:
1343: . snes - the SNES context
1345: Output Parameters:
1346: . sweeps - the number of sweeps of GS to perform.
1348: Level: intermediate
1350: .keywords: SNES, nonlinear, set, Gauss-Siedel
1352: .seealso: SNESSetGS(), SNESGetGS(), SNESSetPC(), SNESSetGSSweeps()
1353: @*/
1354: PetscErrorCode SNESGetGSSweeps(SNES snes, PetscInt * sweeps) {
1356: *sweeps = snes->gssweeps;
1357: return(0);
1358: }
1362: PetscErrorCode SNESPicardComputeFunction(SNES snes,Vec x,Vec f,void *ctx)
1363: {
1366: /* A(x)*x - b(x) */
1367: (*snes->ops->computepjacobian)(snes,x,&snes->jacobian,&snes->jacobian_pre,&snes->matstruct,snes->jacP);
1368: (*snes->ops->computepfunction)(snes,x,f,snes->funP);
1369: VecView(x,PETSC_VIEWER_BINARY_WORLD);
1370: VecView(f,PETSC_VIEWER_BINARY_WORLD);
1371: VecScale(f,-1.0);
1372: MatMultAdd(snes->jacobian_pre,x,f,f);
1373: return(0);
1374: }
1378: PetscErrorCode SNESPicardComputeJacobian(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx)
1379: {
1381: *flag = snes->matstruct;
1382: return(0);
1383: }
1387: /*@C
1388: SNESSetPicard - Use SNES to solve the semilinear-system A(x) x = b(x) via a Picard type iteration (Picard linearization)
1390: Logically Collective on SNES
1392: Input Parameters:
1393: + snes - the SNES context
1394: . r - vector to store function value
1395: . func - function evaluation routine
1396: . jmat - normally the same as mat but you can pass another matrix for which you compute the Jacobian of A(x) x - b(x) (see jmat below)
1397: . mat - matrix to store A
1398: . mfunc - function to compute matrix value
1399: - ctx - [optional] user-defined context for private data for the
1400: function evaluation routine (may be PETSC_NULL)
1402: Calling sequence of func:
1403: $ func (SNES snes,Vec x,Vec f,void *ctx);
1405: + f - function vector
1406: - ctx - optional user-defined function context
1408: Calling sequence of mfunc:
1409: $ mfunc (SNES snes,Vec x,Mat *jmat,Mat *mat,int *flag,void *ctx);
1411: + x - input vector
1412: . jmat - Form Jacobian matrix of A(x) x - b(x) if available, not there is really no reason to use it in this way since then you can just use SNESSetJacobian(),
1413: normally just pass mat in this location
1414: . mat - form A(x) matrix
1415: . flag - flag indicating information about the preconditioner matrix
1416: structure (same as flag in KSPSetOperators()), one of SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER
1417: - ctx - [optional] user-defined Jacobian context
1419: Notes:
1420: One can call SNESSetPicard() or SNESSetFunction() (and possibly SNESSetJacobian()) but cannot call both
1422: $ 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}
1423: $ Note that when an exact solver is used this corresponds to the "classic" Picard A(x^{n}) x^{n+1} = b(x^{n}) iteration.
1425: Run with -snes_mf_operator to solve the system with Newton's method using A(x^{n}) to construct the preconditioner.
1427: We implement the defect correction form of the Picard iteration because it converges much more generally when inexact linear solvers are used then
1428: the direct Picard iteration A(x^n) x^{n+1} = b(x^n)
1430: 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
1431: 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
1432: different please contact us at petsc-dev@mcs.anl.gov and we'll have an entirely new argument :-).
1434: Level: beginner
1436: .keywords: SNES, nonlinear, set, function
1438: .seealso: SNESGetFunction(), SNESSetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESGetPicard(), SNESLineSearchPreCheckPicard()
1439: @*/
1440: PetscErrorCode SNESSetPicard(SNES snes,Vec r,PetscErrorCode (*func)(SNES,Vec,Vec,void*),Mat jmat, Mat mat, PetscErrorCode (*mfunc)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void *ctx)
1441: {
1445: snes->ops->computepfunction = func;
1446: snes->ops->computepjacobian = mfunc;
1447: SNESSetFunction(snes,r,SNESPicardComputeFunction,ctx);
1448: SNESSetJacobian(snes,jmat,mat,SNESPicardComputeJacobian,ctx);
1449: return(0);
1450: }
1454: /*@C
1455: SNESSetComputeInitialGuess - Sets a routine used to compute an initial guess for the problem
1457: Logically Collective on SNES
1459: Input Parameters:
1460: + snes - the SNES context
1461: . func - function evaluation routine
1462: - ctx - [optional] user-defined context for private data for the
1463: function evaluation routine (may be PETSC_NULL)
1465: Calling sequence of func:
1466: $ func (SNES snes,Vec x,void *ctx);
1468: . f - function vector
1469: - ctx - optional user-defined function context
1471: Level: intermediate
1473: .keywords: SNES, nonlinear, set, function
1475: .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian()
1476: @*/
1477: PetscErrorCode SNESSetComputeInitialGuess(SNES snes,PetscErrorCode (*func)(SNES,Vec,void*),void *ctx)
1478: {
1481: if (func) snes->ops->computeinitialguess = func;
1482: if (ctx) snes->initialguessP = ctx;
1483: return(0);
1484: }
1486: /* --------------------------------------------------------------- */
1489: /*@C
1490: SNESGetRhs - Gets the vector for solving F(x) = rhs. If rhs is not set
1491: it assumes a zero right hand side.
1493: Logically Collective on SNES
1495: Input Parameter:
1496: . snes - the SNES context
1498: Output Parameter:
1499: . rhs - the right hand side vector or PETSC_NULL if the right hand side vector is null
1501: Level: intermediate
1503: .keywords: SNES, nonlinear, get, function, right hand side
1505: .seealso: SNESGetSolution(), SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction()
1506: @*/
1507: PetscErrorCode SNESGetRhs(SNES snes,Vec *rhs)
1508: {
1512: *rhs = snes->vec_rhs;
1513: return(0);
1514: }
1518: /*@
1519: SNESComputeFunction - Calls the function that has been set with
1520: SNESSetFunction().
1522: Collective on SNES
1524: Input Parameters:
1525: + snes - the SNES context
1526: - x - input vector
1528: Output Parameter:
1529: . y - function vector, as set by SNESSetFunction()
1531: Notes:
1532: SNESComputeFunction() is typically used within nonlinear solvers
1533: implementations, so most users would not generally call this routine
1534: themselves.
1536: Level: developer
1538: .keywords: SNES, nonlinear, compute, function
1540: .seealso: SNESSetFunction(), SNESGetFunction()
1541: @*/
1542: PetscErrorCode SNESComputeFunction(SNES snes,Vec x,Vec y)
1543: {
1552: VecValidValues(x,2,PETSC_TRUE);
1554: PetscLogEventBegin(SNES_FunctionEval,snes,x,y,0);
1555: if (snes->ops->computefunction) {
1556: PetscStackPush("SNES user function");
1557: (*snes->ops->computefunction)(snes,x,y,snes->funP);
1558: PetscStackPop;
1559: } else if (snes->dm) {
1560: DMComputeFunction(snes->dm,x,y);
1561: } else if (snes->vec_rhs) {
1562: MatMult(snes->jacobian, x, y);
1563: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetFunction() or SNESSetDM() before SNESComputeFunction(), likely called from SNESSolve().");
1564: if (snes->vec_rhs) {
1565: VecAXPY(y,-1.0,snes->vec_rhs);
1566: }
1567: snes->nfuncs++;
1568: PetscLogEventEnd(SNES_FunctionEval,snes,x,y,0);
1569: VecValidValues(y,3,PETSC_FALSE);
1570: return(0);
1571: }
1575: /*@
1576: SNESComputeGS - Calls the Gauss-Seidel function that has been set with
1577: SNESSetGS().
1579: Collective on SNES
1581: Input Parameters:
1582: + snes - the SNES context
1583: . x - input vector
1584: - b - rhs vector
1586: Output Parameter:
1587: . x - new solution vector
1589: Notes:
1590: SNESComputeGS() is typically used within composed nonlinear solver
1591: implementations, so most users would not generally call this routine
1592: themselves.
1594: Level: developer
1596: .keywords: SNES, nonlinear, compute, function
1598: .seealso: SNESSetGS(), SNESComputeFunction()
1599: @*/
1600: PetscErrorCode SNESComputeGS(SNES snes,Vec b,Vec x)
1601: {
1603: PetscInt i;
1611: if (b) VecValidValues(b,2,PETSC_TRUE);
1612: PetscLogEventBegin(SNES_GSEval,snes,x,b,0);
1613: if (snes->ops->computegs) {
1614: for (i = 0; i < snes->gssweeps; i++) {
1615: PetscStackPush("SNES user GS");
1616: (*snes->ops->computegs)(snes,x,b,snes->gsP);
1617: PetscStackPop;
1618: }
1619: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetGS() before SNESComputeGS(), likely called from SNESSolve().");
1620: PetscLogEventEnd(SNES_GSEval,snes,x,b,0);
1621: VecValidValues(x,3,PETSC_FALSE);
1622: return(0);
1623: }
1628: /*@
1629: SNESComputeJacobian - Computes the Jacobian matrix that has been
1630: set with SNESSetJacobian().
1632: Collective on SNES and Mat
1634: Input Parameters:
1635: + snes - the SNES context
1636: - x - input vector
1638: Output Parameters:
1639: + A - Jacobian matrix
1640: . B - optional preconditioning matrix
1641: - flag - flag indicating matrix structure (one of, SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER)
1643: Options Database Keys:
1644: + -snes_lag_preconditioner <lag>
1645: . -snes_lag_jacobian <lag>
1646: . -snes_compare_explicit - Compare the computed Jacobian to the finite difference Jacobian and output the differences
1647: . -snes_compare_explicit_draw - Compare the computed Jacobian to the finite difference Jacobian and draw the result
1648: . -snes_compare_explicit_contour - Compare the computed Jacobian to the finite difference Jacobian and draw a contour plot with the result
1649: . -snes_compare_operator - Make the comparison options above use the operator instead of the preconditioning matrix
1650: . -snes_compare_coloring - Compute the finite differece Jacobian using coloring and display norms of difference
1651: . -snes_compare_coloring_display - Compute the finite differece Jacobian using coloring and display verbose differences
1652: . -snes_compare_coloring_threshold - Display only those matrix entries that differ by more than a given threshold
1653: . -snes_compare_coloring_threshold_atol - Absolute tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold
1654: . -snes_compare_coloring_threshold_rtol - Relative tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold
1655: . -snes_compare_coloring_draw - Compute the finite differece Jacobian using coloring and draw differences
1656: - -snes_compare_coloring_draw_contour - Compute the finite differece Jacobian using coloring and show contours of matrices and differences
1659: Notes:
1660: Most users should not need to explicitly call this routine, as it
1661: is used internally within the nonlinear solvers.
1663: See KSPSetOperators() for important information about setting the
1664: flag parameter.
1666: Level: developer
1668: .keywords: SNES, compute, Jacobian, matrix
1670: .seealso: SNESSetJacobian(), KSPSetOperators(), MatStructure, SNESSetLagPreconditioner(), SNESSetLagJacobian()
1671: @*/
1672: PetscErrorCode SNESComputeJacobian(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *flg)
1673: {
1675: PetscBool flag;
1682: VecValidValues(X,2,PETSC_TRUE);
1683: if (!snes->ops->computejacobian) return(0);
1685: /* make sure that MatAssemblyBegin/End() is called on A matrix if it is matrix free */
1687: if (snes->lagjacobian == -2) {
1688: snes->lagjacobian = -1;
1689: PetscInfo(snes,"Recomputing Jacobian/preconditioner because lag is -2 (means compute Jacobian, but then never again) \n");
1690: } else if (snes->lagjacobian == -1) {
1691: *flg = SAME_PRECONDITIONER;
1692: PetscInfo(snes,"Reusing Jacobian/preconditioner because lag is -1\n");
1693: PetscTypeCompare((PetscObject)*A,MATMFFD,&flag);
1694: if (flag) {
1695: MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
1696: MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
1697: }
1698: return(0);
1699: } else if (snes->lagjacobian > 1 && snes->iter % snes->lagjacobian) {
1700: *flg = SAME_PRECONDITIONER;
1701: PetscInfo2(snes,"Reusing Jacobian/preconditioner because lag is %D and SNES iteration is %D\n",snes->lagjacobian,snes->iter);
1702: PetscTypeCompare((PetscObject)*A,MATMFFD,&flag);
1703: if (flag) {
1704: MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
1705: MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
1706: }
1707: return(0);
1708: }
1710: *flg = DIFFERENT_NONZERO_PATTERN;
1711: PetscLogEventBegin(SNES_JacobianEval,snes,X,*A,*B);
1712: PetscStackPush("SNES user Jacobian function");
1713: (*snes->ops->computejacobian)(snes,X,A,B,flg,snes->jacP);
1714: PetscStackPop;
1715: PetscLogEventEnd(SNES_JacobianEval,snes,X,*A,*B);
1717: if (snes->lagpreconditioner == -2) {
1718: PetscInfo(snes,"Rebuilding preconditioner exactly once since lag is -2\n");
1719: snes->lagpreconditioner = -1;
1720: } else if (snes->lagpreconditioner == -1) {
1721: *flg = SAME_PRECONDITIONER;
1722: PetscInfo(snes,"Reusing preconditioner because lag is -1\n");
1723: } else if (snes->lagpreconditioner > 1 && snes->iter % snes->lagpreconditioner) {
1724: *flg = SAME_PRECONDITIONER;
1725: PetscInfo2(snes,"Reusing preconditioner because lag is %D and SNES iteration is %D\n",snes->lagpreconditioner,snes->iter);
1726: }
1728: /* make sure user returned a correct Jacobian and preconditioner */
1731: {
1732: PetscBool flag = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_operator = PETSC_FALSE;
1733: PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit",&flag,PETSC_NULL);
1734: PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw",&flag_draw,PETSC_NULL);
1735: PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw_contour",&flag_contour,PETSC_NULL);
1736: PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_operator",&flag_operator,PETSC_NULL);
1737: if (flag || flag_draw || flag_contour) {
1738: Mat Bexp_mine = PETSC_NULL,Bexp,FDexp;
1739: MatStructure mstruct;
1740: PetscViewer vdraw,vstdout;
1741: PetscBool flg;
1742: if (flag_operator) {
1743: MatComputeExplicitOperator(*A,&Bexp_mine);
1744: Bexp = Bexp_mine;
1745: } else {
1746: /* See if the preconditioning matrix can be viewed and added directly */
1747: PetscTypeCompareAny((PetscObject)*B,&flg,MATSEQAIJ,MATMPIAIJ,MATSEQDENSE,MATMPIDENSE,MATSEQBAIJ,MATMPIBAIJ,MATSEQSBAIJ,MATMPIBAIJ,"");
1748: if (flg) Bexp = *B;
1749: else {
1750: /* If the "preconditioning" matrix is itself MATSHELL or some other type without direct support */
1751: MatComputeExplicitOperator(*B,&Bexp_mine);
1752: Bexp = Bexp_mine;
1753: }
1754: }
1755: MatConvert(Bexp,MATSAME,MAT_INITIAL_MATRIX,&FDexp);
1756: SNESDefaultComputeJacobian(snes,X,&FDexp,&FDexp,&mstruct,NULL);
1757: PetscViewerASCIIGetStdout(((PetscObject)snes)->comm,&vstdout);
1758: if (flag_draw || flag_contour) {
1759: PetscViewerDrawOpen(((PetscObject)snes)->comm,0,"Explicit Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);
1760: if (flag_contour) {PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);}
1761: } else vdraw = PETSC_NULL;
1762: PetscViewerASCIIPrintf(vstdout,"Explicit %s\n",flag_operator?"Jacobian":"preconditioning Jacobian");
1763: if (flag) {MatView(Bexp,vstdout);}
1764: if (vdraw) {MatView(Bexp,vdraw);}
1765: PetscViewerASCIIPrintf(vstdout,"Finite difference Jacobian\n");
1766: if (flag) {MatView(FDexp,vstdout);}
1767: if (vdraw) {MatView(FDexp,vdraw);}
1768: MatAYPX(FDexp,-1.0,Bexp,SAME_NONZERO_PATTERN);
1769: PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian\n");
1770: if (flag) {MatView(FDexp,vstdout);}
1771: if (vdraw) { /* Always use contour for the difference */
1772: PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);
1773: MatView(FDexp,vdraw);
1774: PetscViewerPopFormat(vdraw);
1775: }
1776: if (flag_contour) {PetscViewerPopFormat(vdraw);}
1777: PetscViewerDestroy(&vdraw);
1778: MatDestroy(&Bexp_mine);
1779: MatDestroy(&FDexp);
1780: }
1781: }
1782: {
1783: PetscBool flag = PETSC_FALSE,flag_display = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_threshold = PETSC_FALSE;
1784: PetscReal threshold_atol = PETSC_SQRT_MACHINE_EPSILON,threshold_rtol = 10*PETSC_SQRT_MACHINE_EPSILON;
1785: PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring",&flag,PETSC_NULL);
1786: PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_display",&flag_display,PETSC_NULL);
1787: PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw",&flag_draw,PETSC_NULL);
1788: PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw_contour",&flag_contour,PETSC_NULL);
1789: PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold",&flag_threshold,PETSC_NULL);
1790: PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_rtol",&threshold_rtol,PETSC_NULL);
1791: PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_atol",&threshold_atol,PETSC_NULL);
1792: if (flag || flag_display || flag_draw || flag_contour || flag_threshold) {
1793: Mat Bfd;
1794: MatStructure mstruct;
1795: PetscViewer vdraw,vstdout;
1796: ISColoring iscoloring;
1797: MatFDColoring matfdcoloring;
1798: PetscErrorCode (*func)(SNES,Vec,Vec,void*);
1799: void *funcctx;
1800: PetscReal norm1,norm2,normmax;
1802: MatDuplicate(*B,MAT_DO_NOT_COPY_VALUES,&Bfd);
1803: MatGetColoring(Bfd,MATCOLORINGSL,&iscoloring);
1804: MatFDColoringCreate(Bfd,iscoloring,&matfdcoloring);
1805: ISColoringDestroy(&iscoloring);
1807: /* This method of getting the function is currently unreliable since it doesn't work for DM local functions. */
1808: SNESGetFunction(snes,PETSC_NULL,&func,&funcctx);
1809: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode(*)(void))func,funcctx);
1810: PetscObjectSetOptionsPrefix((PetscObject)matfdcoloring,((PetscObject)snes)->prefix);
1811: PetscObjectAppendOptionsPrefix((PetscObject)matfdcoloring,"coloring_");
1812: MatFDColoringSetFromOptions(matfdcoloring);
1813: MatFDColoringApply(Bfd,matfdcoloring,X,&mstruct,snes);
1814: MatFDColoringDestroy(&matfdcoloring);
1816: PetscViewerASCIIGetStdout(((PetscObject)snes)->comm,&vstdout);
1817: if (flag_draw || flag_contour) {
1818: PetscViewerDrawOpen(((PetscObject)snes)->comm,0,"Colored Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);
1819: if (flag_contour) {PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);}
1820: } else vdraw = PETSC_NULL;
1821: PetscViewerASCIIPrintf(vstdout,"Explicit preconditioning Jacobian\n");
1822: if (flag_display) {MatView(*B,vstdout);}
1823: if (vdraw) {MatView(*B,vdraw);}
1824: PetscViewerASCIIPrintf(vstdout,"Colored Finite difference Jacobian\n");
1825: if (flag_display) {MatView(Bfd,vstdout);}
1826: if (vdraw) {MatView(Bfd,vdraw);}
1827: MatAYPX(Bfd,-1.0,*B,SAME_NONZERO_PATTERN);
1828: MatNorm(Bfd,NORM_1,&norm1);
1829: MatNorm(Bfd,NORM_FROBENIUS,&norm2);
1830: MatNorm(Bfd,NORM_MAX,&normmax);
1831: PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian, norm1=%G normFrob=%G normmax=%G\n",norm1,norm2,normmax);
1832: if (flag_display) {MatView(Bfd,vstdout);}
1833: if (vdraw) { /* Always use contour for the difference */
1834: PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);
1835: MatView(Bfd,vdraw);
1836: PetscViewerPopFormat(vdraw);
1837: }
1838: if (flag_contour) {PetscViewerPopFormat(vdraw);}
1840: if (flag_threshold) {
1841: PetscInt bs,rstart,rend,i;
1842: MatGetBlockSize(*B,&bs);
1843: MatGetOwnershipRange(*B,&rstart,&rend);
1844: for (i=rstart; i<rend; i++) {
1845: const PetscScalar *ba,*ca;
1846: const PetscInt *bj,*cj;
1847: PetscInt bn,cn,j,maxentrycol = -1,maxdiffcol = -1,maxrdiffcol = -1;
1848: PetscReal maxentry = 0,maxdiff = 0,maxrdiff = 0;
1849: MatGetRow(*B,i,&bn,&bj,&ba);
1850: MatGetRow(Bfd,i,&cn,&cj,&ca);
1851: if (bn != cn) SETERRQ(((PetscObject)*A)->comm,PETSC_ERR_PLIB,"Unexpected different nonzero pattern in -snes_compare_coloring_threshold");
1852: for (j=0; j<bn; j++) {
1853: PetscReal rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j]));
1854: if (PetscAbsScalar(ba[j]) > PetscAbs(maxentry)) {
1855: maxentrycol = bj[j];
1856: maxentry = PetscRealPart(ba[j]);
1857: }
1858: if (PetscAbsScalar(ca[j]) > PetscAbs(maxdiff)) {
1859: maxdiffcol = bj[j];
1860: maxdiff = PetscRealPart(ca[j]);
1861: }
1862: if (rdiff > maxrdiff) {
1863: maxrdiffcol = bj[j];
1864: maxrdiff = rdiff;
1865: }
1866: }
1867: if (maxrdiff > 1) {
1868: PetscViewerASCIIPrintf(vstdout,"row %D (maxentry=%G at %D, maxdiff=%G at %D, maxrdiff=%G at %D):",i,maxentry,maxentrycol,maxdiff,maxdiffcol,maxrdiff,maxrdiffcol);
1869: for (j=0; j<bn; j++) {
1870: PetscReal rdiff;
1871: rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j]));
1872: if (rdiff > 1) {
1873: PetscViewerASCIIPrintf(vstdout," (%D,%G:%G)",bj[j],PetscRealPart(ba[j]),PetscRealPart(ca[j]));
1874: }
1875: }
1876: PetscViewerASCIIPrintf(vstdout,"\n",i,maxentry,maxdiff,maxrdiff);
1877: }
1878: MatRestoreRow(*B,i,&bn,&bj,&ba);
1879: MatRestoreRow(Bfd,i,&cn,&cj,&ca);
1880: }
1881: }
1882: PetscViewerDestroy(&vdraw);
1883: MatDestroy(&Bfd);
1884: }
1885: }
1886: return(0);
1887: }
1891: /*@C
1892: SNESSetJacobian - Sets the function to compute Jacobian as well as the
1893: location to store the matrix.
1895: Logically Collective on SNES and Mat
1897: Input Parameters:
1898: + snes - the SNES context
1899: . A - Jacobian matrix
1900: . B - preconditioner matrix (usually same as the Jacobian)
1901: . func - Jacobian evaluation routine (if PETSC_NULL then SNES retains any previously set value)
1902: - ctx - [optional] user-defined context for private data for the
1903: Jacobian evaluation routine (may be PETSC_NULL) (if PETSC_NULL then SNES retains any previously set value)
1905: Calling sequence of func:
1906: $ func (SNES snes,Vec x,Mat *A,Mat *B,int *flag,void *ctx);
1908: + x - input vector
1909: . A - Jacobian matrix
1910: . B - preconditioner matrix, usually the same as A
1911: . flag - flag indicating information about the preconditioner matrix
1912: structure (same as flag in KSPSetOperators()), one of SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER
1913: - ctx - [optional] user-defined Jacobian context
1915: Notes:
1916: See KSPSetOperators() for important information about setting the flag
1917: output parameter in the routine func(). Be sure to read this information!
1919: The routine func() takes Mat * as the matrix arguments rather than Mat.
1920: This allows the Jacobian evaluation routine to replace A and/or B with a
1921: completely new new matrix structure (not just different matrix elements)
1922: when appropriate, for instance, if the nonzero structure is changing
1923: throughout the global iterations.
1925: If the A matrix and B matrix are different you must call MatAssemblyBegin/End() on
1926: each matrix.
1928: If using SNESDefaultComputeJacobianColor() to assemble a Jacobian, the ctx argument
1929: must be a MatFDColoring.
1931: Other defect-correction schemes can be used by computing a different matrix in place of the Jacobian. One common
1932: example is to use the "Picard linearization" which only differentiates through the highest order parts of each term.
1934: Level: beginner
1936: .keywords: SNES, nonlinear, set, Jacobian, matrix
1938: .seealso: KSPSetOperators(), SNESSetFunction(), MatMFFDComputeJacobian(), SNESDefaultComputeJacobianColor(), MatStructure
1939: @*/
1940: PetscErrorCode SNESSetJacobian(SNES snes,Mat A,Mat B,PetscErrorCode (*func)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void *ctx)
1941: {
1950: if (func) snes->ops->computejacobian = func;
1951: if (ctx) snes->jacP = ctx;
1952: if (A) {
1953: PetscObjectReference((PetscObject)A);
1954: MatDestroy(&snes->jacobian);
1955: snes->jacobian = A;
1956: }
1957: if (B) {
1958: PetscObjectReference((PetscObject)B);
1959: MatDestroy(&snes->jacobian_pre);
1960: snes->jacobian_pre = B;
1961: }
1962: return(0);
1963: }
1967: /*@C
1968: SNESGetJacobian - Returns the Jacobian matrix and optionally the user
1969: provided context for evaluating the Jacobian.
1971: Not Collective, but Mat object will be parallel if SNES object is
1973: Input Parameter:
1974: . snes - the nonlinear solver context
1976: Output Parameters:
1977: + A - location to stash Jacobian matrix (or PETSC_NULL)
1978: . B - location to stash preconditioner matrix (or PETSC_NULL)
1979: . func - location to put Jacobian function (or PETSC_NULL)
1980: - ctx - location to stash Jacobian ctx (or PETSC_NULL)
1982: Level: advanced
1984: .seealso: SNESSetJacobian(), SNESComputeJacobian()
1985: @*/
1986: PetscErrorCode SNESGetJacobian(SNES snes,Mat *A,Mat *B,PetscErrorCode (**func)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void **ctx)
1987: {
1990: if (A) *A = snes->jacobian;
1991: if (B) *B = snes->jacobian_pre;
1992: if (func) *func = snes->ops->computejacobian;
1993: if (ctx) *ctx = snes->jacP;
1994: return(0);
1995: }
1997: /* ----- Routines to initialize and destroy a nonlinear solver ---- */
2001: /*@
2002: SNESSetUp - Sets up the internal data structures for the later use
2003: of a nonlinear solver.
2005: Collective on SNES
2007: Input Parameters:
2008: . snes - the SNES context
2010: Notes:
2011: For basic use of the SNES solvers the user need not explicitly call
2012: SNESSetUp(), since these actions will automatically occur during
2013: the call to SNESSolve(). However, if one wishes to control this
2014: phase separately, SNESSetUp() should be called after SNESCreate()
2015: and optional routines of the form SNESSetXXX(), but before SNESSolve().
2017: Level: advanced
2019: .keywords: SNES, nonlinear, setup
2021: .seealso: SNESCreate(), SNESSolve(), SNESDestroy()
2022: @*/
2023: PetscErrorCode SNESSetUp(SNES snes)
2024: {
2029: if (snes->setupcalled) return(0);
2031: if (!((PetscObject)snes)->type_name) {
2032: SNESSetType(snes,SNESLS);
2033: }
2035: SNESGetFunction(snes,&snes->vec_func,PETSC_NULL,PETSC_NULL);
2036: if (snes->vec_func == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be function vector");
2037: if (snes->vec_rhs == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be right hand side vector");
2039: if (!snes->vec_sol_update /* && snes->vec_sol */) {
2040: VecDuplicate(snes->vec_sol,&snes->vec_sol_update);
2041: PetscLogObjectParent(snes,snes->vec_sol_update);
2042: }
2044: if (!snes->ops->computejacobian && snes->dm) {
2045: Mat J,B;
2046: DMCreateMatrix(snes->dm,MATAIJ,&B);
2047: if (snes->mf_operator) {
2048: MatCreateSNESMF(snes,&J);
2049: MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);
2050: MatSetFromOptions(J);
2051: } else {
2052: J = B;
2053: PetscObjectReference((PetscObject)J);
2054: }
2055: SNESSetJacobian(snes,J,B,SNESDMComputeJacobian,PETSC_NULL);
2056: MatDestroy(&J);
2057: MatDestroy(&B);
2058: } else if (!snes->jacobian && snes->ops->computejacobian == MatMFFDComputeJacobian) {
2059: Mat J;
2060: MatCreateSNESMF(snes,&J);
2061: MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);
2062: MatSetFromOptions(J);
2063: SNESSetJacobian(snes,J,J,MatMFFDComputeJacobian,snes->funP);
2064: MatDestroy(&J);
2065: } else if (snes->dm && snes->mf_operator && !snes->jacobian_pre && !snes->jacobian) {
2066: Mat J,B;
2067: MatCreateSNESMF(snes,&J);
2068: MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);
2069: MatSetFromOptions(J);
2070: DMCreateMatrix(snes->dm,MATAIJ,&B);
2071: SNESSetJacobian(snes,J,B,SNESDMComputeJacobian,snes->funP);
2072: MatDestroy(&J);
2073: MatDestroy(&B);
2074: } else if (snes->dm && !snes->jacobian_pre){
2075: Mat J;
2076: DMCreateMatrix(snes->dm,MATAIJ,&J);
2077: SNESSetJacobian(snes,J,J,PETSC_NULL,PETSC_NULL);
2078: MatDestroy(&J);
2079: }
2080: if (!snes->ops->computefunction && !snes->dm) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must provide a residual function with SNESSetFunction() or call SNESSetDM()");
2081: if (!snes->vec_func) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must provide a vector when calling SNESSetFunction() or call SNESSetDM()");
2083: if (!snes->ksp) {SNESGetKSP(snes, &snes->ksp);}
2085: if (snes->ops->usercompute && !snes->user) {
2086: (*snes->ops->usercompute)(snes,(void**)&snes->user);
2087: }
2089: if (snes->ops->setup) {
2090: (*snes->ops->setup)(snes);
2091: }
2093: snes->setupcalled = PETSC_TRUE;
2094: return(0);
2095: }
2099: /*@
2100: SNESReset - Resets a SNES context to the snessetupcalled = 0 state and removes any allocated Vecs and Mats
2102: Collective on SNES
2104: Input Parameter:
2105: . snes - iterative context obtained from SNESCreate()
2107: Level: intermediate
2109: Notes: Also calls the application context destroy routine set with SNESSetComputeApplicationContext()
2111: .keywords: SNES, destroy
2113: .seealso: SNESCreate(), SNESSetUp(), SNESSolve()
2114: @*/
2115: PetscErrorCode SNESReset(SNES snes)
2116: {
2121: if (snes->ops->userdestroy && snes->user) {
2122: (*snes->ops->userdestroy)((void**)&snes->user);
2123: snes->user = PETSC_NULL;
2124: }
2125: if (snes->pc) {
2126: SNESReset(snes->pc);
2127: }
2129: if (snes->ops->reset) {
2130: (*snes->ops->reset)(snes);
2131: }
2132: if (snes->ksp) {KSPReset(snes->ksp);}
2133: VecDestroy(&snes->vec_rhs);
2134: VecDestroy(&snes->vec_sol);
2135: VecDestroy(&snes->vec_sol_update);
2136: VecDestroy(&snes->vec_func);
2137: MatDestroy(&snes->jacobian);
2138: MatDestroy(&snes->jacobian_pre);
2139: VecDestroyVecs(snes->nwork,&snes->work);
2140: VecDestroyVecs(snes->nvwork,&snes->vwork);
2141: snes->nwork = snes->nvwork = 0;
2142: snes->setupcalled = PETSC_FALSE;
2143: return(0);
2144: }
2148: /*@
2149: SNESDestroy - Destroys the nonlinear solver context that was created
2150: with SNESCreate().
2152: Collective on SNES
2154: Input Parameter:
2155: . snes - the SNES context
2157: Level: beginner
2159: .keywords: SNES, nonlinear, destroy
2161: .seealso: SNESCreate(), SNESSolve()
2162: @*/
2163: PetscErrorCode SNESDestroy(SNES *snes)
2164: {
2168: if (!*snes) return(0);
2170: if (--((PetscObject)(*snes))->refct > 0) {*snes = 0; return(0);}
2172: SNESReset((*snes));
2173: SNESDestroy(&(*snes)->pc);
2175: /* if memory was published with AMS then destroy it */
2176: PetscObjectDepublish((*snes));
2177: if ((*snes)->ops->destroy) {(*((*snes))->ops->destroy)((*snes));}
2179: DMDestroy(&(*snes)->dm);
2180: KSPDestroy(&(*snes)->ksp);
2182: PetscFree((*snes)->kspconvctx);
2183: if ((*snes)->ops->convergeddestroy) {
2184: (*(*snes)->ops->convergeddestroy)((*snes)->cnvP);
2185: }
2186: if ((*snes)->conv_malloc) {
2187: PetscFree((*snes)->conv_hist);
2188: PetscFree((*snes)->conv_hist_its);
2189: }
2190: PetscViewerDestroy(&(*snes)->ls_monitor);
2191: SNESMonitorCancel((*snes));
2192: PetscHeaderDestroy(snes);
2193: return(0);
2194: }
2196: /* ----------- Routines to set solver parameters ---------- */
2200: /*@
2201: SNESSetLagPreconditioner - Determines when the preconditioner is rebuilt in the nonlinear solve.
2203: Logically Collective on SNES
2205: Input Parameters:
2206: + snes - the SNES context
2207: - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time
2208: the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that
2210: Options Database Keys:
2211: . -snes_lag_preconditioner <lag>
2213: Notes:
2214: The default is 1
2215: The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1
2216: If -1 is used before the very first nonlinear solve the preconditioner is still built because there is no previous preconditioner to use
2218: Level: intermediate
2220: .keywords: SNES, nonlinear, set, convergence, tolerances
2222: .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian()
2224: @*/
2225: PetscErrorCode SNESSetLagPreconditioner(SNES snes,PetscInt lag)
2226: {
2229: if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater");
2230: if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0");
2232: snes->lagpreconditioner = lag;
2233: return(0);
2234: }
2238: /*@
2239: SNESSetGridSequence - sets the number of steps of grid sequencing that SNES does
2241: Logically Collective on SNES
2243: Input Parameters:
2244: + snes - the SNES context
2245: - steps - the number of refinements to do, defaults to 0
2247: Options Database Keys:
2248: . -snes_grid_sequence <steps>
2250: Level: intermediate
2252: Notes:
2253: Use SNESGetSolution() to extract the fine grid solution after grid sequencing.
2255: .keywords: SNES, nonlinear, set, convergence, tolerances
2257: .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian()
2259: @*/
2260: PetscErrorCode SNESSetGridSequence(SNES snes,PetscInt steps)
2261: {
2265: snes->gridsequence = steps;
2266: return(0);
2267: }
2271: /*@
2272: SNESGetLagPreconditioner - Indicates how often the preconditioner is rebuilt
2274: Not Collective
2276: Input Parameter:
2277: . snes - the SNES context
2278:
2279: Output Parameter:
2280: . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time
2281: the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that
2283: Options Database Keys:
2284: . -snes_lag_preconditioner <lag>
2286: Notes:
2287: The default is 1
2288: The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1
2290: Level: intermediate
2292: .keywords: SNES, nonlinear, set, convergence, tolerances
2294: .seealso: SNESSetTrustRegionTolerance(), SNESSetLagPreconditioner()
2296: @*/
2297: PetscErrorCode SNESGetLagPreconditioner(SNES snes,PetscInt *lag)
2298: {
2301: *lag = snes->lagpreconditioner;
2302: return(0);
2303: }
2307: /*@
2308: SNESSetLagJacobian - Determines when the Jacobian is rebuilt in the nonlinear solve. See SNESSetLagPreconditioner() for determining how
2309: often the preconditioner is rebuilt.
2311: Logically Collective on SNES
2313: Input Parameters:
2314: + snes - the SNES context
2315: - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time
2316: the Jacobian is built etc. -2 means rebuild at next chance but then never again
2318: Options Database Keys:
2319: . -snes_lag_jacobian <lag>
2321: Notes:
2322: The default is 1
2323: The Jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1
2324: 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
2325: at the next Newton step but never again (unless it is reset to another value)
2327: Level: intermediate
2329: .keywords: SNES, nonlinear, set, convergence, tolerances
2331: .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagPreconditioner(), SNESGetLagJacobian()
2333: @*/
2334: PetscErrorCode SNESSetLagJacobian(SNES snes,PetscInt lag)
2335: {
2338: if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater");
2339: if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0");
2341: snes->lagjacobian = lag;
2342: return(0);
2343: }
2347: /*@
2348: SNESGetLagJacobian - Indicates how often the Jacobian is rebuilt. See SNESGetLagPreconditioner() to determine when the preconditioner is rebuilt
2350: Not Collective
2352: Input Parameter:
2353: . snes - the SNES context
2354:
2355: Output Parameter:
2356: . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time
2357: the Jacobian is built etc.
2359: Options Database Keys:
2360: . -snes_lag_jacobian <lag>
2362: Notes:
2363: The default is 1
2364: The jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1
2366: Level: intermediate
2368: .keywords: SNES, nonlinear, set, convergence, tolerances
2370: .seealso: SNESSetTrustRegionTolerance(), SNESSetLagJacobian(), SNESSetLagPreconditioner(), SNESGetLagPreconditioner()
2372: @*/
2373: PetscErrorCode SNESGetLagJacobian(SNES snes,PetscInt *lag)
2374: {
2377: *lag = snes->lagjacobian;
2378: return(0);
2379: }
2383: /*@
2384: SNESSetTolerances - Sets various parameters used in convergence tests.
2386: Logically Collective on SNES
2388: Input Parameters:
2389: + snes - the SNES context
2390: . abstol - absolute convergence tolerance
2391: . rtol - relative convergence tolerance
2392: . stol - convergence tolerance in terms of the norm
2393: of the change in the solution between steps
2394: . maxit - maximum number of iterations
2395: - maxf - maximum number of function evaluations
2397: Options Database Keys:
2398: + -snes_atol <abstol> - Sets abstol
2399: . -snes_rtol <rtol> - Sets rtol
2400: . -snes_stol <stol> - Sets stol
2401: . -snes_max_it <maxit> - Sets maxit
2402: - -snes_max_funcs <maxf> - Sets maxf
2404: Notes:
2405: The default maximum number of iterations is 50.
2406: The default maximum number of function evaluations is 1000.
2408: Level: intermediate
2410: .keywords: SNES, nonlinear, set, convergence, tolerances
2412: .seealso: SNESSetTrustRegionTolerance()
2413: @*/
2414: PetscErrorCode SNESSetTolerances(SNES snes,PetscReal abstol,PetscReal rtol,PetscReal stol,PetscInt maxit,PetscInt maxf)
2415: {
2424: if (abstol != PETSC_DEFAULT) {
2425: if (abstol < 0.0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %G must be non-negative",abstol);
2426: snes->abstol = abstol;
2427: }
2428: if (rtol != PETSC_DEFAULT) {
2429: if (rtol < 0.0 || 1.0 <= rtol) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Relative tolerance %G must be non-negative and less than 1.0",rtol);
2430: snes->rtol = rtol;
2431: }
2432: if (stol != PETSC_DEFAULT) {
2433: if (stol < 0.0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Step tolerance %G must be non-negative",stol);
2434: snes->xtol = stol;
2435: }
2436: if (maxit != PETSC_DEFAULT) {
2437: if (maxit < 0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxit);
2438: snes->max_its = maxit;
2439: }
2440: if (maxf != PETSC_DEFAULT) {
2441: if (maxf < 0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of function evaluations %D must be non-negative",maxf);
2442: snes->max_funcs = maxf;
2443: }
2444: return(0);
2445: }
2449: /*@
2450: SNESGetTolerances - Gets various parameters used in convergence tests.
2452: Not Collective
2454: Input Parameters:
2455: + snes - the SNES context
2456: . atol - absolute convergence tolerance
2457: . rtol - relative convergence tolerance
2458: . stol - convergence tolerance in terms of the norm
2459: of the change in the solution between steps
2460: . maxit - maximum number of iterations
2461: - maxf - maximum number of function evaluations
2463: Notes:
2464: The user can specify PETSC_NULL for any parameter that is not needed.
2466: Level: intermediate
2468: .keywords: SNES, nonlinear, get, convergence, tolerances
2470: .seealso: SNESSetTolerances()
2471: @*/
2472: PetscErrorCode SNESGetTolerances(SNES snes,PetscReal *atol,PetscReal *rtol,PetscReal *stol,PetscInt *maxit,PetscInt *maxf)
2473: {
2476: if (atol) *atol = snes->abstol;
2477: if (rtol) *rtol = snes->rtol;
2478: if (stol) *stol = snes->xtol;
2479: if (maxit) *maxit = snes->max_its;
2480: if (maxf) *maxf = snes->max_funcs;
2481: return(0);
2482: }
2486: /*@
2487: SNESSetTrustRegionTolerance - Sets the trust region parameter tolerance.
2489: Logically Collective on SNES
2491: Input Parameters:
2492: + snes - the SNES context
2493: - tol - tolerance
2494:
2495: Options Database Key:
2496: . -snes_trtol <tol> - Sets tol
2498: Level: intermediate
2500: .keywords: SNES, nonlinear, set, trust region, tolerance
2502: .seealso: SNESSetTolerances()
2503: @*/
2504: PetscErrorCode SNESSetTrustRegionTolerance(SNES snes,PetscReal tol)
2505: {
2509: snes->deltatol = tol;
2510: return(0);
2511: }
2513: /*
2514: Duplicate the lg monitors for SNES from KSP; for some reason with
2515: dynamic libraries things don't work under Sun4 if we just use
2516: macros instead of functions
2517: */
2520: PetscErrorCode SNESMonitorLG(SNES snes,PetscInt it,PetscReal norm,void *ctx)
2521: {
2526: KSPMonitorLG((KSP)snes,it,norm,ctx);
2527: return(0);
2528: }
2532: PetscErrorCode SNESMonitorLGCreate(const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *draw)
2533: {
2537: KSPMonitorLGCreate(host,label,x,y,m,n,draw);
2538: return(0);
2539: }
2543: PetscErrorCode SNESMonitorLGDestroy(PetscDrawLG *draw)
2544: {
2548: KSPMonitorLGDestroy(draw);
2549: return(0);
2550: }
2552: extern PetscErrorCode SNESMonitorRange_Private(SNES,PetscInt,PetscReal*);
2555: PetscErrorCode SNESMonitorLGRange(SNES snes,PetscInt n,PetscReal rnorm,void *monctx)
2556: {
2557: PetscDrawLG lg;
2558: PetscErrorCode ierr;
2559: PetscReal x,y,per;
2560: PetscViewer v = (PetscViewer)monctx;
2561: static PetscReal prev; /* should be in the context */
2562: PetscDraw draw;
2564: if (!monctx) {
2565: MPI_Comm comm;
2567: PetscObjectGetComm((PetscObject)snes,&comm);
2568: v = PETSC_VIEWER_DRAW_(comm);
2569: }
2570: PetscViewerDrawGetDrawLG(v,0,&lg);
2571: if (!n) {PetscDrawLGReset(lg);}
2572: PetscDrawLGGetDraw(lg,&draw);
2573: PetscDrawSetTitle(draw,"Residual norm");
2574: x = (PetscReal) n;
2575: if (rnorm > 0.0) y = log10(rnorm); else y = -15.0;
2576: PetscDrawLGAddPoint(lg,&x,&y);
2577: if (n < 20 || !(n % 5)) {
2578: PetscDrawLGDraw(lg);
2579: }
2581: PetscViewerDrawGetDrawLG(v,1,&lg);
2582: if (!n) {PetscDrawLGReset(lg);}
2583: PetscDrawLGGetDraw(lg,&draw);
2584: PetscDrawSetTitle(draw,"% elemts > .2*max elemt");
2585: SNESMonitorRange_Private(snes,n,&per);
2586: x = (PetscReal) n;
2587: y = 100.0*per;
2588: PetscDrawLGAddPoint(lg,&x,&y);
2589: if (n < 20 || !(n % 5)) {
2590: PetscDrawLGDraw(lg);
2591: }
2593: PetscViewerDrawGetDrawLG(v,2,&lg);
2594: if (!n) {prev = rnorm;PetscDrawLGReset(lg);}
2595: PetscDrawLGGetDraw(lg,&draw);
2596: PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm");
2597: x = (PetscReal) n;
2598: y = (prev - rnorm)/prev;
2599: PetscDrawLGAddPoint(lg,&x,&y);
2600: if (n < 20 || !(n % 5)) {
2601: PetscDrawLGDraw(lg);
2602: }
2604: PetscViewerDrawGetDrawLG(v,3,&lg);
2605: if (!n) {PetscDrawLGReset(lg);}
2606: PetscDrawLGGetDraw(lg,&draw);
2607: PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm*(% > .2 max)");
2608: x = (PetscReal) n;
2609: y = (prev - rnorm)/(prev*per);
2610: if (n > 2) { /*skip initial crazy value */
2611: PetscDrawLGAddPoint(lg,&x,&y);
2612: }
2613: if (n < 20 || !(n % 5)) {
2614: PetscDrawLGDraw(lg);
2615: }
2616: prev = rnorm;
2617: return(0);
2618: }
2622: PetscErrorCode SNESMonitorLGRangeCreate(const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *draw)
2623: {
2627: KSPMonitorLGCreate(host,label,x,y,m,n,draw);
2628: return(0);
2629: }
2633: PetscErrorCode SNESMonitorLGRangeDestroy(PetscDrawLG *draw)
2634: {
2638: KSPMonitorLGDestroy(draw);
2639: return(0);
2640: }
2644: /*@
2645: SNESMonitor - runs the user provided monitor routines, if they exist
2647: Collective on SNES
2649: Input Parameters:
2650: + snes - nonlinear solver context obtained from SNESCreate()
2651: . iter - iteration number
2652: - rnorm - relative norm of the residual
2654: Notes:
2655: This routine is called by the SNES implementations.
2656: It does not typically need to be called by the user.
2658: Level: developer
2660: .seealso: SNESMonitorSet()
2661: @*/
2662: PetscErrorCode SNESMonitor(SNES snes,PetscInt iter,PetscReal rnorm)
2663: {
2665: PetscInt i,n = snes->numbermonitors;
2668: for (i=0; i<n; i++) {
2669: (*snes->monitor[i])(snes,iter,rnorm,snes->monitorcontext[i]);
2670: }
2671: return(0);
2672: }
2674: /* ------------ Routines to set performance monitoring options ----------- */
2678: /*@C
2679: SNESMonitorSet - Sets an ADDITIONAL function that is to be used at every
2680: iteration of the nonlinear solver to display the iteration's
2681: progress.
2683: Logically Collective on SNES
2685: Input Parameters:
2686: + snes - the SNES context
2687: . func - monitoring routine
2688: . mctx - [optional] user-defined context for private data for the
2689: monitor routine (use PETSC_NULL if no context is desired)
2690: - monitordestroy - [optional] routine that frees monitor context
2691: (may be PETSC_NULL)
2693: Calling sequence of func:
2694: $ int func(SNES snes,PetscInt its, PetscReal norm,void *mctx)
2696: + snes - the SNES context
2697: . its - iteration number
2698: . norm - 2-norm function value (may be estimated)
2699: - mctx - [optional] monitoring context
2701: Options Database Keys:
2702: + -snes_monitor - sets SNESMonitorDefault()
2703: . -snes_monitor_draw - sets line graph monitor,
2704: uses SNESMonitorLGCreate()
2705: - -snes_monitor_cancel - cancels all monitors that have
2706: been hardwired into a code by
2707: calls to SNESMonitorSet(), but
2708: does not cancel those set via
2709: the options database.
2711: Notes:
2712: Several different monitoring routines may be set by calling
2713: SNESMonitorSet() multiple times; all will be called in the
2714: order in which they were set.
2716: Fortran notes: Only a single monitor function can be set for each SNES object
2718: Level: intermediate
2720: .keywords: SNES, nonlinear, set, monitor
2722: .seealso: SNESMonitorDefault(), SNESMonitorCancel()
2723: @*/
2724: PetscErrorCode SNESMonitorSet(SNES snes,PetscErrorCode (*monitor)(SNES,PetscInt,PetscReal,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**))
2725: {
2726: PetscInt i;
2731: if (snes->numbermonitors >= MAXSNESMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set");
2732: for (i=0; i<snes->numbermonitors;i++) {
2733: if (monitor == snes->monitor[i] && monitordestroy == snes->monitordestroy[i] && mctx == snes->monitorcontext[i]) {
2734: if (monitordestroy) {
2735: (*monitordestroy)(&mctx);
2736: }
2737: return(0);
2738: }
2739: }
2740: snes->monitor[snes->numbermonitors] = monitor;
2741: snes->monitordestroy[snes->numbermonitors] = monitordestroy;
2742: snes->monitorcontext[snes->numbermonitors++] = (void*)mctx;
2743: return(0);
2744: }
2748: /*@C
2749: SNESMonitorCancel - Clears all the monitor functions for a SNES object.
2751: Logically Collective on SNES
2753: Input Parameters:
2754: . snes - the SNES context
2756: Options Database Key:
2757: . -snes_monitor_cancel - cancels all monitors that have been hardwired
2758: into a code by calls to SNESMonitorSet(), but does not cancel those
2759: set via the options database
2761: Notes:
2762: There is no way to clear one specific monitor from a SNES object.
2764: Level: intermediate
2766: .keywords: SNES, nonlinear, set, monitor
2768: .seealso: SNESMonitorDefault(), SNESMonitorSet()
2769: @*/
2770: PetscErrorCode SNESMonitorCancel(SNES snes)
2771: {
2773: PetscInt i;
2777: for (i=0; i<snes->numbermonitors; i++) {
2778: if (snes->monitordestroy[i]) {
2779: (*snes->monitordestroy[i])(&snes->monitorcontext[i]);
2780: }
2781: }
2782: snes->numbermonitors = 0;
2783: return(0);
2784: }
2788: /*@C
2789: SNESSetConvergenceTest - Sets the function that is to be used
2790: to test for convergence of the nonlinear iterative solution.
2792: Logically Collective on SNES
2794: Input Parameters:
2795: + snes - the SNES context
2796: . func - routine to test for convergence
2797: . cctx - [optional] context for private data for the convergence routine (may be PETSC_NULL)
2798: - destroy - [optional] destructor for the context (may be PETSC_NULL; PETSC_NULL_FUNCTION in Fortran)
2800: Calling sequence of func:
2801: $ PetscErrorCode func (SNES snes,PetscInt it,PetscReal xnorm,PetscReal gnorm,PetscReal f,SNESConvergedReason *reason,void *cctx)
2803: + snes - the SNES context
2804: . it - current iteration (0 is the first and is before any Newton step)
2805: . cctx - [optional] convergence context
2806: . reason - reason for convergence/divergence
2807: . xnorm - 2-norm of current iterate
2808: . gnorm - 2-norm of current step
2809: - f - 2-norm of function
2811: Level: advanced
2813: .keywords: SNES, nonlinear, set, convergence, test
2815: .seealso: SNESDefaultConverged(), SNESSkipConverged()
2816: @*/
2817: PetscErrorCode SNESSetConvergenceTest(SNES snes,PetscErrorCode (*func)(SNES,PetscInt,PetscReal,PetscReal,PetscReal,SNESConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*))
2818: {
2823: if (!func) func = SNESSkipConverged;
2824: if (snes->ops->convergeddestroy) {
2825: (*snes->ops->convergeddestroy)(snes->cnvP);
2826: }
2827: snes->ops->converged = func;
2828: snes->ops->convergeddestroy = destroy;
2829: snes->cnvP = cctx;
2830: return(0);
2831: }
2835: /*@
2836: SNESGetConvergedReason - Gets the reason the SNES iteration was stopped.
2838: Not Collective
2840: Input Parameter:
2841: . snes - the SNES context
2843: Output Parameter:
2844: . reason - negative value indicates diverged, positive value converged, see SNESConvergedReason or the
2845: manual pages for the individual convergence tests for complete lists
2847: Level: intermediate
2849: Notes: Can only be called after the call the SNESSolve() is complete.
2851: .keywords: SNES, nonlinear, set, convergence, test
2853: .seealso: SNESSetConvergenceTest(), SNESConvergedReason
2854: @*/
2855: PetscErrorCode SNESGetConvergedReason(SNES snes,SNESConvergedReason *reason)
2856: {
2860: *reason = snes->reason;
2861: return(0);
2862: }
2866: /*@
2867: SNESSetConvergenceHistory - Sets the array used to hold the convergence history.
2869: Logically Collective on SNES
2871: Input Parameters:
2872: + snes - iterative context obtained from SNESCreate()
2873: . a - array to hold history, this array will contain the function norms computed at each step
2874: . its - integer array holds the number of linear iterations for each solve.
2875: . na - size of a and its
2876: - reset - PETSC_TRUE indicates each new nonlinear solve resets the history counter to zero,
2877: else it continues storing new values for new nonlinear solves after the old ones
2879: Notes:
2880: If 'a' and 'its' are PETSC_NULL then space is allocated for the history. If 'na' PETSC_DECIDE or PETSC_DEFAULT then a
2881: default array of length 10000 is allocated.
2883: This routine is useful, e.g., when running a code for purposes
2884: of accurate performance monitoring, when no I/O should be done
2885: during the section of code that is being timed.
2887: Level: intermediate
2889: .keywords: SNES, set, convergence, history
2891: .seealso: SNESGetConvergenceHistory()
2893: @*/
2894: PetscErrorCode SNESSetConvergenceHistory(SNES snes,PetscReal a[],PetscInt its[],PetscInt na,PetscBool reset)
2895: {
2902: if (na == PETSC_DECIDE || na == PETSC_DEFAULT || !a) {
2903: if (na == PETSC_DECIDE || na == PETSC_DEFAULT) na = 1000;
2904: PetscMalloc(na*sizeof(PetscReal),&a);
2905: PetscMalloc(na*sizeof(PetscInt),&its);
2906: snes->conv_malloc = PETSC_TRUE;
2907: }
2908: snes->conv_hist = a;
2909: snes->conv_hist_its = its;
2910: snes->conv_hist_max = na;
2911: snes->conv_hist_len = 0;
2912: snes->conv_hist_reset = reset;
2913: return(0);
2914: }
2916: #if defined(PETSC_HAVE_MATLAB_ENGINE)
2917: #include <engine.h> /* MATLAB include file */
2918: #include <mex.h> /* MATLAB include file */
2919: EXTERN_C_BEGIN
2922: mxArray *SNESGetConvergenceHistoryMatlab(SNES snes)
2923: {
2924: mxArray *mat;
2925: PetscInt i;
2926: PetscReal *ar;
2929: mat = mxCreateDoubleMatrix(snes->conv_hist_len,1,mxREAL);
2930: ar = (PetscReal*) mxGetData(mat);
2931: for (i=0; i<snes->conv_hist_len; i++) {
2932: ar[i] = snes->conv_hist[i];
2933: }
2934: PetscFunctionReturn(mat);
2935: }
2936: EXTERN_C_END
2937: #endif
2942: /*@C
2943: SNESGetConvergenceHistory - Gets the array used to hold the convergence history.
2945: Not Collective
2947: Input Parameter:
2948: . snes - iterative context obtained from SNESCreate()
2950: Output Parameters:
2951: . a - array to hold history
2952: . its - integer array holds the number of linear iterations (or
2953: negative if not converged) for each solve.
2954: - na - size of a and its
2956: Notes:
2957: The calling sequence for this routine in Fortran is
2958: $ call SNESGetConvergenceHistory(SNES snes, integer na, integer ierr)
2960: This routine is useful, e.g., when running a code for purposes
2961: of accurate performance monitoring, when no I/O should be done
2962: during the section of code that is being timed.
2964: Level: intermediate
2966: .keywords: SNES, get, convergence, history
2968: .seealso: SNESSetConvergencHistory()
2970: @*/
2971: PetscErrorCode SNESGetConvergenceHistory(SNES snes,PetscReal *a[],PetscInt *its[],PetscInt *na)
2972: {
2975: if (a) *a = snes->conv_hist;
2976: if (its) *its = snes->conv_hist_its;
2977: if (na) *na = snes->conv_hist_len;
2978: return(0);
2979: }
2983: /*@C
2984: SNESSetUpdate - Sets the general-purpose update function called
2985: at the beginning of every iteration of the nonlinear solve. Specifically
2986: it is called just before the Jacobian is "evaluated".
2988: Logically Collective on SNES
2990: Input Parameters:
2991: . snes - The nonlinear solver context
2992: . func - The function
2994: Calling sequence of func:
2995: . func (SNES snes, PetscInt step);
2997: . step - The current step of the iteration
2999: Level: advanced
3001: 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()
3002: This is not used by most users.
3004: .keywords: SNES, update
3006: .seealso SNESDefaultUpdate(), SNESSetJacobian(), SNESSolve()
3007: @*/
3008: PetscErrorCode SNESSetUpdate(SNES snes, PetscErrorCode (*func)(SNES, PetscInt))
3009: {
3012: snes->ops->update = func;
3013: return(0);
3014: }
3018: /*@
3019: SNESDefaultUpdate - The default update function which does nothing.
3021: Not collective
3023: Input Parameters:
3024: . snes - The nonlinear solver context
3025: . step - The current step of the iteration
3027: Level: intermediate
3029: .keywords: SNES, update
3030: .seealso SNESSetUpdate(), SNESDefaultRhsBC(), SNESDefaultShortolutionBC()
3031: @*/
3032: PetscErrorCode SNESDefaultUpdate(SNES snes, PetscInt step)
3033: {
3035: return(0);
3036: }
3040: /*
3041: SNESScaleStep_Private - Scales a step so that its length is less than the
3042: positive parameter delta.
3044: Input Parameters:
3045: + snes - the SNES context
3046: . y - approximate solution of linear system
3047: . fnorm - 2-norm of current function
3048: - delta - trust region size
3050: Output Parameters:
3051: + gpnorm - predicted function norm at the new point, assuming local
3052: linearization. The value is zero if the step lies within the trust
3053: region, and exceeds zero otherwise.
3054: - ynorm - 2-norm of the step
3056: Note:
3057: For non-trust region methods such as SNESLS, the parameter delta
3058: is set to be the maximum allowable step size.
3060: .keywords: SNES, nonlinear, scale, step
3061: */
3062: PetscErrorCode SNESScaleStep_Private(SNES snes,Vec y,PetscReal *fnorm,PetscReal *delta,PetscReal *gpnorm,PetscReal *ynorm)
3063: {
3064: PetscReal nrm;
3065: PetscScalar cnorm;
3073: VecNorm(y,NORM_2,&nrm);
3074: if (nrm > *delta) {
3075: nrm = *delta/nrm;
3076: *gpnorm = (1.0 - nrm)*(*fnorm);
3077: cnorm = nrm;
3078: VecScale(y,cnorm);
3079: *ynorm = *delta;
3080: } else {
3081: *gpnorm = 0.0;
3082: *ynorm = nrm;
3083: }
3084: return(0);
3085: }
3089: /*@C
3090: SNESSolve - Solves a nonlinear system F(x) = b.
3091: Call SNESSolve() after calling SNESCreate() and optional routines of the form SNESSetXXX().
3093: Collective on SNES
3095: Input Parameters:
3096: + snes - the SNES context
3097: . b - the constant part of the equation F(x) = b, or PETSC_NULL to use zero.
3098: - x - the solution vector.
3100: Notes:
3101: The user should initialize the vector,x, with the initial guess
3102: for the nonlinear solve prior to calling SNESSolve. In particular,
3103: to employ an initial guess of zero, the user should explicitly set
3104: this vector to zero by calling VecSet().
3106: Level: beginner
3108: .keywords: SNES, nonlinear, solve
3110: .seealso: SNESCreate(), SNESDestroy(), SNESSetFunction(), SNESSetJacobian(), SNESSetGridSequence(), SNESGetSolution()
3111: @*/
3112: PetscErrorCode SNESSolve(SNES snes,Vec b,Vec x)
3113: {
3115: PetscBool flg;
3116: char filename[PETSC_MAX_PATH_LEN];
3117: PetscViewer viewer;
3118: PetscInt grid;
3119: Vec xcreated = PETSC_NULL;
3128: if (!x && snes->dm) {
3129: DMCreateGlobalVector(snes->dm,&xcreated);
3130: x = xcreated;
3131: }
3133: for (grid=0; grid<snes->gridsequence; grid++) {PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm));}
3134: for (grid=0; grid<snes->gridsequence+1; grid++) {
3136: /* set solution vector */
3137: if (!grid) {PetscObjectReference((PetscObject)x);}
3138: VecDestroy(&snes->vec_sol);
3139: snes->vec_sol = x;
3140: /* set afine vector if provided */
3141: if (b) { PetscObjectReference((PetscObject)b); }
3142: VecDestroy(&snes->vec_rhs);
3143: snes->vec_rhs = b;
3145: SNESSetUp(snes);
3147: if (!grid) {
3148: if (snes->ops->computeinitialguess) {
3149: (*snes->ops->computeinitialguess)(snes,snes->vec_sol,snes->initialguessP);
3150: } else if (snes->dm) {
3151: PetscBool ig;
3152: DMHasInitialGuess(snes->dm,&ig);
3153: if (ig) {
3154: DMComputeInitialGuess(snes->dm,snes->vec_sol);
3155: }
3156: }
3157: }
3159: if (snes->conv_hist_reset) snes->conv_hist_len = 0;
3160: snes->nfuncs = 0; snes->linear_its = 0; snes->numFailures = 0;
3162: PetscLogEventBegin(SNES_Solve,snes,0,0,0);
3163: (*snes->ops->solve)(snes);
3164: PetscLogEventEnd(SNES_Solve,snes,0,0,0);
3165: if (snes->domainerror){
3166: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
3167: snes->domainerror = PETSC_FALSE;
3168: }
3169: if (!snes->reason) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
3171: PetscOptionsGetString(((PetscObject)snes)->prefix,"-snes_view",filename,PETSC_MAX_PATH_LEN,&flg);
3172: if (flg && !PetscPreLoadingOn) {
3173: PetscViewerASCIIOpen(((PetscObject)snes)->comm,filename,&viewer);
3174: SNESView(snes,viewer);
3175: PetscViewerDestroy(&viewer);
3176: }
3177:
3178: flg = PETSC_FALSE;
3179: PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_test_local_min",&flg,PETSC_NULL);
3180: if (flg && !PetscPreLoadingOn) { SNESTestLocalMin(snes); }
3181: if (snes->printreason) {
3182: PetscViewerASCIIAddTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),((PetscObject)snes)->tablevel);
3183: if (snes->reason > 0) {
3184: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),"Nonlinear solve converged due to %s\n",SNESConvergedReasons[snes->reason]);
3185: } else {
3186: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),"Nonlinear solve did not converge due to %s\n",SNESConvergedReasons[snes->reason]);
3187: }
3188: PetscViewerASCIISubtractTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),((PetscObject)snes)->tablevel);
3189: }
3191: flg = PETSC_FALSE;
3192: PetscOptionsGetString(((PetscObject)snes)->prefix,"-snes_view_solution_vtk",filename,PETSC_MAX_PATH_LEN,&flg);
3193: if (flg) {
3194: PetscViewer viewer;
3195: PetscViewerCreate(((PetscObject)snes)->comm,&viewer);
3196: PetscViewerSetType(viewer,PETSCVIEWERVTK);
3197: PetscViewerFileSetName(viewer,filename);
3198: VecView(snes->vec_sol,viewer);
3199: PetscViewerDestroy(&viewer);
3200: }
3202: if (snes->errorifnotconverged && snes->reason < 0) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged");
3203: if (grid < snes->gridsequence) {
3204: DM fine;
3205: Vec xnew;
3206: Mat interp;
3208: DMRefine(snes->dm,((PetscObject)snes)->comm,&fine);
3209: DMCreateInterpolation(snes->dm,fine,&interp,PETSC_NULL);
3210: DMCreateGlobalVector(fine,&xnew);
3211: MatInterpolate(interp,x,xnew);
3212: MatDestroy(&interp);
3213: x = xnew;
3215: SNESReset(snes);
3216: SNESSetDM(snes,fine);
3217: DMDestroy(&fine);
3218: PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm));
3219: }
3220: }
3221: VecDestroy(&xcreated);
3222: return(0);
3223: }
3225: /* --------- Internal routines for SNES Package --------- */
3229: /*@C
3230: SNESSetType - Sets the method for the nonlinear solver.
3232: Collective on SNES
3234: Input Parameters:
3235: + snes - the SNES context
3236: - type - a known method
3238: Options Database Key:
3239: . -snes_type <type> - Sets the method; use -help for a list
3240: of available methods (for instance, ls or tr)
3242: Notes:
3243: See "petsc/include/petscsnes.h" for available methods (for instance)
3244: + SNESLS - Newton's method with line search
3245: (systems of nonlinear equations)
3246: . SNESTR - Newton's method with trust region
3247: (systems of nonlinear equations)
3249: Normally, it is best to use the SNESSetFromOptions() command and then
3250: set the SNES solver type from the options database rather than by using
3251: this routine. Using the options database provides the user with
3252: maximum flexibility in evaluating the many nonlinear solvers.
3253: The SNESSetType() routine is provided for those situations where it
3254: is necessary to set the nonlinear solver independently of the command
3255: line or options database. This might be the case, for example, when
3256: the choice of solver changes during the execution of the program,
3257: and the user's application is taking responsibility for choosing the
3258: appropriate method.
3260: Level: intermediate
3262: .keywords: SNES, set, type
3264: .seealso: SNESType, SNESCreate()
3266: @*/
3267: PetscErrorCode SNESSetType(SNES snes,const SNESType type)
3268: {
3269: PetscErrorCode ierr,(*r)(SNES);
3270: PetscBool match;
3276: PetscTypeCompare((PetscObject)snes,type,&match);
3277: if (match) return(0);
3279: PetscFListFind(SNESList,((PetscObject)snes)->comm,type,PETSC_TRUE,(void (**)(void)) &r);
3280: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested SNES type %s",type);
3281: /* Destroy the previous private SNES context */
3282: if (snes->ops->destroy) {
3283: (*(snes)->ops->destroy)(snes);
3284: snes->ops->destroy = PETSC_NULL;
3285: }
3286: /* Reinitialize function pointers in SNESOps structure */
3287: snes->ops->setup = 0;
3288: snes->ops->solve = 0;
3289: snes->ops->view = 0;
3290: snes->ops->setfromoptions = 0;
3291: snes->ops->destroy = 0;
3292: /* Call the SNESCreate_XXX routine for this particular Nonlinear solver */
3293: snes->setupcalled = PETSC_FALSE;
3294: PetscObjectChangeTypeName((PetscObject)snes,type);
3295: (*r)(snes);
3296: #if defined(PETSC_HAVE_AMS)
3297: if (PetscAMSPublishAll) {
3298: PetscObjectAMSPublish((PetscObject)snes);
3299: }
3300: #endif
3301: return(0);
3302: }
3305: /* --------------------------------------------------------------------- */
3308: /*@
3309: SNESRegisterDestroy - Frees the list of nonlinear solvers that were
3310: registered by SNESRegisterDynamic().
3312: Not Collective
3314: Level: advanced
3316: .keywords: SNES, nonlinear, register, destroy
3318: .seealso: SNESRegisterAll(), SNESRegisterAll()
3319: @*/
3320: PetscErrorCode SNESRegisterDestroy(void)
3321: {
3325: PetscFListDestroy(&SNESList);
3326: SNESRegisterAllCalled = PETSC_FALSE;
3327: return(0);
3328: }
3332: /*@C
3333: SNESGetType - Gets the SNES method type and name (as a string).
3335: Not Collective
3337: Input Parameter:
3338: . snes - nonlinear solver context
3340: Output Parameter:
3341: . type - SNES method (a character string)
3343: Level: intermediate
3345: .keywords: SNES, nonlinear, get, type, name
3346: @*/
3347: PetscErrorCode SNESGetType(SNES snes,const SNESType *type)
3348: {
3352: *type = ((PetscObject)snes)->type_name;
3353: return(0);
3354: }
3358: /*@
3359: SNESGetSolution - Returns the vector where the approximate solution is
3360: stored. This is the fine grid solution when using SNESSetGridSequence().
3362: Not Collective, but Vec is parallel if SNES is parallel
3364: Input Parameter:
3365: . snes - the SNES context
3367: Output Parameter:
3368: . x - the solution
3370: Level: intermediate
3372: .keywords: SNES, nonlinear, get, solution
3374: .seealso: SNESGetSolutionUpdate(), SNESGetFunction()
3375: @*/
3376: PetscErrorCode SNESGetSolution(SNES snes,Vec *x)
3377: {
3381: *x = snes->vec_sol;
3382: return(0);
3383: }
3387: /*@
3388: SNESGetSolutionUpdate - Returns the vector where the solution update is
3389: stored.
3391: Not Collective, but Vec is parallel if SNES is parallel
3393: Input Parameter:
3394: . snes - the SNES context
3396: Output Parameter:
3397: . x - the solution update
3399: Level: advanced
3401: .keywords: SNES, nonlinear, get, solution, update
3403: .seealso: SNESGetSolution(), SNESGetFunction()
3404: @*/
3405: PetscErrorCode SNESGetSolutionUpdate(SNES snes,Vec *x)
3406: {
3410: *x = snes->vec_sol_update;
3411: return(0);
3412: }
3416: /*@C
3417: SNESGetFunction - Returns the vector where the function is stored.
3419: Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet.
3421: Input Parameter:
3422: . snes - the SNES context
3424: Output Parameter:
3425: + r - the function (or PETSC_NULL)
3426: . func - the function (or PETSC_NULL)
3427: - ctx - the function context (or PETSC_NULL)
3429: Level: advanced
3431: .keywords: SNES, nonlinear, get, function
3433: .seealso: SNESSetFunction(), SNESGetSolution()
3434: @*/
3435: PetscErrorCode SNESGetFunction(SNES snes,Vec *r,PetscErrorCode (**func)(SNES,Vec,Vec,void*),void **ctx)
3436: {
3441: if (r) {
3442: if (!snes->vec_func) {
3443: if (snes->vec_rhs) {
3444: VecDuplicate(snes->vec_rhs,&snes->vec_func);
3445: } else if (snes->vec_sol) {
3446: VecDuplicate(snes->vec_sol,&snes->vec_func);
3447: } else if (snes->dm) {
3448: DMCreateGlobalVector(snes->dm,&snes->vec_func);
3449: }
3450: }
3451: *r = snes->vec_func;
3452: }
3453: if (func) *func = snes->ops->computefunction;
3454: if (ctx) *ctx = snes->funP;
3455: return(0);
3456: }
3458: /*@C
3459: SNESGetGS - Returns the GS function and context.
3461: Input Parameter:
3462: . snes - the SNES context
3464: Output Parameter:
3465: + gsfunc - the function (or PETSC_NULL)
3466: - ctx - the function context (or PETSC_NULL)
3468: Level: advanced
3470: .keywords: SNES, nonlinear, get, function
3472: .seealso: SNESSetGS(), SNESGetFunction()
3473: @*/
3477: PetscErrorCode SNESGetGS (SNES snes, PetscErrorCode(**func)(SNES, Vec, Vec, void*), void ** ctx)
3478: {
3481: if (func) *func = snes->ops->computegs;
3482: if (ctx) *ctx = snes->funP;
3483: return(0);
3484: }
3488: /*@C
3489: SNESSetOptionsPrefix - Sets the prefix used for searching for all
3490: SNES options in the database.
3492: Logically Collective on SNES
3494: Input Parameter:
3495: + snes - the SNES context
3496: - prefix - the prefix to prepend to all option names
3498: Notes:
3499: A hyphen (-) must NOT be given at the beginning of the prefix name.
3500: The first character of all runtime options is AUTOMATICALLY the hyphen.
3502: Level: advanced
3504: .keywords: SNES, set, options, prefix, database
3506: .seealso: SNESSetFromOptions()
3507: @*/
3508: PetscErrorCode SNESSetOptionsPrefix(SNES snes,const char prefix[])
3509: {
3514: PetscObjectSetOptionsPrefix((PetscObject)snes,prefix);
3515: if (!snes->ksp) {SNESGetKSP(snes,&snes->ksp);}
3516: KSPSetOptionsPrefix(snes->ksp,prefix);
3517: return(0);
3518: }
3522: /*@C
3523: SNESAppendOptionsPrefix - Appends to the prefix used for searching for all
3524: SNES options in the database.
3526: Logically Collective on SNES
3528: Input Parameters:
3529: + snes - the SNES context
3530: - prefix - the prefix to prepend to all option names
3532: Notes:
3533: A hyphen (-) must NOT be given at the beginning of the prefix name.
3534: The first character of all runtime options is AUTOMATICALLY the hyphen.
3536: Level: advanced
3538: .keywords: SNES, append, options, prefix, database
3540: .seealso: SNESGetOptionsPrefix()
3541: @*/
3542: PetscErrorCode SNESAppendOptionsPrefix(SNES snes,const char prefix[])
3543: {
3545:
3548: PetscObjectAppendOptionsPrefix((PetscObject)snes,prefix);
3549: if (!snes->ksp) {SNESGetKSP(snes,&snes->ksp);}
3550: KSPAppendOptionsPrefix(snes->ksp,prefix);
3551: return(0);
3552: }
3556: /*@C
3557: SNESGetOptionsPrefix - Sets the prefix used for searching for all
3558: SNES options in the database.
3560: Not Collective
3562: Input Parameter:
3563: . snes - the SNES context
3565: Output Parameter:
3566: . prefix - pointer to the prefix string used
3568: Notes: On the fortran side, the user should pass in a string 'prefix' of
3569: sufficient length to hold the prefix.
3571: Level: advanced
3573: .keywords: SNES, get, options, prefix, database
3575: .seealso: SNESAppendOptionsPrefix()
3576: @*/
3577: PetscErrorCode SNESGetOptionsPrefix(SNES snes,const char *prefix[])
3578: {
3583: PetscObjectGetOptionsPrefix((PetscObject)snes,prefix);
3584: return(0);
3585: }
3590: /*@C
3591: SNESRegister - See SNESRegisterDynamic()
3593: Level: advanced
3594: @*/
3595: PetscErrorCode SNESRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(SNES))
3596: {
3597: char fullname[PETSC_MAX_PATH_LEN];
3601: PetscFListConcat(path,name,fullname);
3602: PetscFListAdd(&SNESList,sname,fullname,(void (*)(void))function);
3603: return(0);
3604: }
3608: PetscErrorCode SNESTestLocalMin(SNES snes)
3609: {
3611: PetscInt N,i,j;
3612: Vec u,uh,fh;
3613: PetscScalar value;
3614: PetscReal norm;
3617: SNESGetSolution(snes,&u);
3618: VecDuplicate(u,&uh);
3619: VecDuplicate(u,&fh);
3621: /* currently only works for sequential */
3622: PetscPrintf(PETSC_COMM_WORLD,"Testing FormFunction() for local min\n");
3623: VecGetSize(u,&N);
3624: for (i=0; i<N; i++) {
3625: VecCopy(u,uh);
3626: PetscPrintf(PETSC_COMM_WORLD,"i = %D\n",i);
3627: for (j=-10; j<11; j++) {
3628: value = PetscSign(j)*exp(PetscAbs(j)-10.0);
3629: VecSetValue(uh,i,value,ADD_VALUES);
3630: SNESComputeFunction(snes,uh,fh);
3631: VecNorm(fh,NORM_2,&norm);
3632: PetscPrintf(PETSC_COMM_WORLD," j norm %D %18.16e\n",j,norm);
3633: value = -value;
3634: VecSetValue(uh,i,value,ADD_VALUES);
3635: }
3636: }
3637: VecDestroy(&uh);
3638: VecDestroy(&fh);
3639: return(0);
3640: }
3644: /*@
3645: SNESKSPSetUseEW - Sets SNES use Eisenstat-Walker method for
3646: computing relative tolerance for linear solvers within an inexact
3647: Newton method.
3649: Logically Collective on SNES
3651: Input Parameters:
3652: + snes - SNES context
3653: - flag - PETSC_TRUE or PETSC_FALSE
3655: Options Database:
3656: + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence
3657: . -snes_ksp_ew_version ver - version of Eisenstat-Walker method
3658: . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0
3659: . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax
3660: . -snes_ksp_ew_gamma <gamma> - Sets gamma
3661: . -snes_ksp_ew_alpha <alpha> - Sets alpha
3662: . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2
3663: - -snes_ksp_ew_threshold <threshold> - Sets threshold
3665: Notes:
3666: Currently, the default is to use a constant relative tolerance for
3667: the inner linear solvers. Alternatively, one can use the
3668: Eisenstat-Walker method, where the relative convergence tolerance
3669: is reset at each Newton iteration according progress of the nonlinear
3670: solver.
3672: Level: advanced
3674: Reference:
3675: S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an
3676: inexact Newton method", SISC 17 (1), pp.16-32, 1996.
3678: .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton
3680: .seealso: SNESKSPGetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW()
3681: @*/
3682: PetscErrorCode SNESKSPSetUseEW(SNES snes,PetscBool flag)
3683: {
3687: snes->ksp_ewconv = flag;
3688: return(0);
3689: }
3693: /*@
3694: SNESKSPGetUseEW - Gets if SNES is using Eisenstat-Walker method
3695: for computing relative tolerance for linear solvers within an
3696: inexact Newton method.
3698: Not Collective
3700: Input Parameter:
3701: . snes - SNES context
3703: Output Parameter:
3704: . flag - PETSC_TRUE or PETSC_FALSE
3706: Notes:
3707: Currently, the default is to use a constant relative tolerance for
3708: the inner linear solvers. Alternatively, one can use the
3709: Eisenstat-Walker method, where the relative convergence tolerance
3710: is reset at each Newton iteration according progress of the nonlinear
3711: solver.
3713: Level: advanced
3715: Reference:
3716: S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an
3717: inexact Newton method", SISC 17 (1), pp.16-32, 1996.
3719: .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton
3721: .seealso: SNESKSPSetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW()
3722: @*/
3723: PetscErrorCode SNESKSPGetUseEW(SNES snes, PetscBool *flag)
3724: {
3728: *flag = snes->ksp_ewconv;
3729: return(0);
3730: }
3734: /*@
3735: SNESKSPSetParametersEW - Sets parameters for Eisenstat-Walker
3736: convergence criteria for the linear solvers within an inexact
3737: Newton method.
3739: Logically Collective on SNES
3740:
3741: Input Parameters:
3742: + snes - SNES context
3743: . version - version 1, 2 (default is 2) or 3
3744: . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1)
3745: . rtol_max - maximum relative tolerance (0 <= rtol_max < 1)
3746: . gamma - multiplicative factor for version 2 rtol computation
3747: (0 <= gamma2 <= 1)
3748: . alpha - power for version 2 rtol computation (1 < alpha <= 2)
3749: . alpha2 - power for safeguard
3750: - threshold - threshold for imposing safeguard (0 < threshold < 1)
3752: Note:
3753: Version 3 was contributed by Luis Chacon, June 2006.
3755: Use PETSC_DEFAULT to retain the default for any of the parameters.
3757: Level: advanced
3759: Reference:
3760: S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an
3761: inexact Newton method", Utah State University Math. Stat. Dept. Res.
3762: Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput.
3764: .keywords: SNES, KSP, Eisenstat, Walker, set, parameters
3766: .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPGetParametersEW()
3767: @*/
3768: PetscErrorCode SNESKSPSetParametersEW(SNES snes,PetscInt version,PetscReal rtol_0,PetscReal rtol_max,
3769: PetscReal gamma,PetscReal alpha,PetscReal alpha2,PetscReal threshold)
3770: {
3771: SNESKSPEW *kctx;
3774: kctx = (SNESKSPEW*)snes->kspconvctx;
3775: if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing");
3784: if (version != PETSC_DEFAULT) kctx->version = version;
3785: if (rtol_0 != PETSC_DEFAULT) kctx->rtol_0 = rtol_0;
3786: if (rtol_max != PETSC_DEFAULT) kctx->rtol_max = rtol_max;
3787: if (gamma != PETSC_DEFAULT) kctx->gamma = gamma;
3788: if (alpha != PETSC_DEFAULT) kctx->alpha = alpha;
3789: if (alpha2 != PETSC_DEFAULT) kctx->alpha2 = alpha2;
3790: if (threshold != PETSC_DEFAULT) kctx->threshold = threshold;
3791:
3792: if (kctx->version < 1 || kctx->version > 3) {
3793: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 and 3 are supported: %D",kctx->version);
3794: }
3795: if (kctx->rtol_0 < 0.0 || kctx->rtol_0 >= 1.0) {
3796: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= rtol_0 < 1.0: %G",kctx->rtol_0);
3797: }
3798: if (kctx->rtol_max < 0.0 || kctx->rtol_max >= 1.0) {
3799: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= rtol_max (%G) < 1.0\n",kctx->rtol_max);
3800: }
3801: if (kctx->gamma < 0.0 || kctx->gamma > 1.0) {
3802: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= gamma (%G) <= 1.0\n",kctx->gamma);
3803: }
3804: if (kctx->alpha <= 1.0 || kctx->alpha > 2.0) {
3805: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"1.0 < alpha (%G) <= 2.0\n",kctx->alpha);
3806: }
3807: if (kctx->threshold <= 0.0 || kctx->threshold >= 1.0) {
3808: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 < threshold (%G) < 1.0\n",kctx->threshold);
3809: }
3810: return(0);
3811: }
3815: /*@
3816: SNESKSPGetParametersEW - Gets parameters for Eisenstat-Walker
3817: convergence criteria for the linear solvers within an inexact
3818: Newton method.
3820: Not Collective
3821:
3822: Input Parameters:
3823: snes - SNES context
3825: Output Parameters:
3826: + version - version 1, 2 (default is 2) or 3
3827: . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1)
3828: . rtol_max - maximum relative tolerance (0 <= rtol_max < 1)
3829: . gamma - multiplicative factor for version 2 rtol computation
3830: (0 <= gamma2 <= 1)
3831: . alpha - power for version 2 rtol computation (1 < alpha <= 2)
3832: . alpha2 - power for safeguard
3833: - threshold - threshold for imposing safeguard (0 < threshold < 1)
3835: Level: advanced
3837: .keywords: SNES, KSP, Eisenstat, Walker, get, parameters
3839: .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPSetParametersEW()
3840: @*/
3841: PetscErrorCode SNESKSPGetParametersEW(SNES snes,PetscInt *version,PetscReal *rtol_0,PetscReal *rtol_max,
3842: PetscReal *gamma,PetscReal *alpha,PetscReal *alpha2,PetscReal *threshold)
3843: {
3844: SNESKSPEW *kctx;
3847: kctx = (SNESKSPEW*)snes->kspconvctx;
3848: if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing");
3849: if(version) *version = kctx->version;
3850: if(rtol_0) *rtol_0 = kctx->rtol_0;
3851: if(rtol_max) *rtol_max = kctx->rtol_max;
3852: if(gamma) *gamma = kctx->gamma;
3853: if(alpha) *alpha = kctx->alpha;
3854: if(alpha2) *alpha2 = kctx->alpha2;
3855: if(threshold) *threshold = kctx->threshold;
3856: return(0);
3857: }
3861: static PetscErrorCode SNESKSPEW_PreSolve(SNES snes, KSP ksp, Vec b, Vec x)
3862: {
3864: SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx;
3865: PetscReal rtol=PETSC_DEFAULT,stol;
3868: if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context exists");
3869: if (!snes->iter) { /* first time in, so use the original user rtol */
3870: rtol = kctx->rtol_0;
3871: } else {
3872: if (kctx->version == 1) {
3873: rtol = (snes->norm - kctx->lresid_last)/kctx->norm_last;
3874: if (rtol < 0.0) rtol = -rtol;
3875: stol = pow(kctx->rtol_last,kctx->alpha2);
3876: if (stol > kctx->threshold) rtol = PetscMax(rtol,stol);
3877: } else if (kctx->version == 2) {
3878: rtol = kctx->gamma * pow(snes->norm/kctx->norm_last,kctx->alpha);
3879: stol = kctx->gamma * pow(kctx->rtol_last,kctx->alpha);
3880: if (stol > kctx->threshold) rtol = PetscMax(rtol,stol);
3881: } else if (kctx->version == 3) {/* contributed by Luis Chacon, June 2006. */
3882: rtol = kctx->gamma * pow(snes->norm/kctx->norm_last,kctx->alpha);
3883: /* safeguard: avoid sharp decrease of rtol */
3884: stol = kctx->gamma*pow(kctx->rtol_last,kctx->alpha);
3885: stol = PetscMax(rtol,stol);
3886: rtol = PetscMin(kctx->rtol_0,stol);
3887: /* safeguard: avoid oversolving */
3888: stol = kctx->gamma*(snes->ttol)/snes->norm;
3889: stol = PetscMax(rtol,stol);
3890: rtol = PetscMin(kctx->rtol_0,stol);
3891: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 or 3 are supported: %D",kctx->version);
3892: }
3893: /* safeguard: avoid rtol greater than one */
3894: rtol = PetscMin(rtol,kctx->rtol_max);
3895: KSPSetTolerances(ksp,rtol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
3896: PetscInfo3(snes,"iter %D, Eisenstat-Walker (version %D) KSP rtol=%G\n",snes->iter,kctx->version,rtol);
3897: return(0);
3898: }
3902: static PetscErrorCode SNESKSPEW_PostSolve(SNES snes, KSP ksp, Vec b, Vec x)
3903: {
3905: SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx;
3906: PCSide pcside;
3907: Vec lres;
3910: if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context exists");
3911: KSPGetTolerances(ksp,&kctx->rtol_last,0,0,0);
3912: SNESGetFunctionNorm(snes,&kctx->norm_last);
3913: if (kctx->version == 1) {
3914: KSPGetPCSide(ksp,&pcside);
3915: if (pcside == PC_RIGHT) { /* XXX Should we also test KSP_UNPRECONDITIONED_NORM ? */
3916: /* KSP residual is true linear residual */
3917: KSPGetResidualNorm(ksp,&kctx->lresid_last);
3918: } else {
3919: /* KSP residual is preconditioned residual */
3920: /* compute true linear residual norm */
3921: VecDuplicate(b,&lres);
3922: MatMult(snes->jacobian,x,lres);
3923: VecAYPX(lres,-1.0,b);
3924: VecNorm(lres,NORM_2,&kctx->lresid_last);
3925: VecDestroy(&lres);
3926: }
3927: }
3928: return(0);
3929: }
3933: PetscErrorCode SNES_KSPSolve(SNES snes, KSP ksp, Vec b, Vec x)
3934: {
3938: if (snes->ksp_ewconv) { SNESKSPEW_PreSolve(snes,ksp,b,x); }
3939: KSPSolve(ksp,b,x);
3940: if (snes->ksp_ewconv) { SNESKSPEW_PostSolve(snes,ksp,b,x); }
3941: return(0);
3942: }
3946: /*@
3947: SNESSetDM - Sets the DM that may be used by some preconditioners
3949: Logically Collective on SNES
3951: Input Parameters:
3952: + snes - the preconditioner context
3953: - dm - the dm
3955: Level: intermediate
3958: .seealso: SNESGetDM(), KSPSetDM(), KSPGetDM()
3959: @*/
3960: PetscErrorCode SNESSetDM(SNES snes,DM dm)
3961: {
3963: KSP ksp;
3967: if (dm) {PetscObjectReference((PetscObject)dm);}
3968: DMDestroy(&snes->dm);
3969: snes->dm = dm;
3970: SNESGetKSP(snes,&ksp);
3971: KSPSetDM(ksp,dm);
3972: KSPSetDMActive(ksp,PETSC_FALSE);
3973: if (snes->pc) {
3974: SNESSetDM(snes->pc, snes->dm);
3975: }
3976: return(0);
3977: }
3981: /*@
3982: SNESGetDM - Gets the DM that may be used by some preconditioners
3984: Not Collective but DM obtained is parallel on SNES
3986: Input Parameter:
3987: . snes - the preconditioner context
3989: Output Parameter:
3990: . dm - the dm
3992: Level: intermediate
3995: .seealso: SNESSetDM(), KSPSetDM(), KSPGetDM()
3996: @*/
3997: PetscErrorCode SNESGetDM(SNES snes,DM *dm)
3998: {
4001: *dm = snes->dm;
4002: return(0);
4003: }
4007: /*@
4008: SNESSetPC - Sets the nonlinear preconditioner to be used.
4010: Collective on SNES
4012: Input Parameters:
4013: + snes - iterative context obtained from SNESCreate()
4014: - pc - the preconditioner object
4016: Notes:
4017: Use SNESGetPC() to retrieve the preconditioner context (for example,
4018: to configure it using the API).
4020: Level: developer
4022: .keywords: SNES, set, precondition
4023: .seealso: SNESGetPC()
4024: @*/
4025: PetscErrorCode SNESSetPC(SNES snes, SNES pc)
4026: {
4033: PetscObjectReference((PetscObject) pc);
4034: SNESDestroy(&snes->pc);
4035: snes->pc = pc;
4036: PetscLogObjectParent(snes, snes->pc);
4037: return(0);
4038: }
4042: /*@
4043: SNESGetPC - Returns a pointer to the nonlinear preconditioning context set with SNESSetPC().
4045: Not Collective
4047: Input Parameter:
4048: . snes - iterative context obtained from SNESCreate()
4050: Output Parameter:
4051: . pc - preconditioner context
4053: Level: developer
4055: .keywords: SNES, get, preconditioner
4056: .seealso: SNESSetPC()
4057: @*/
4058: PetscErrorCode SNESGetPC(SNES snes, SNES *pc)
4059: {
4065: if (!snes->pc) {
4066: SNESCreate(((PetscObject) snes)->comm, &snes->pc);
4067: PetscObjectIncrementTabLevel((PetscObject) snes->pc, (PetscObject) snes, 1);
4068: PetscLogObjectParent(snes, snes->pc);
4069: }
4070: *pc = snes->pc;
4071: return(0);
4072: }
4074: #if defined(PETSC_HAVE_MATLAB_ENGINE)
4075: #include <mex.h>
4077: typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext;
4081: /*
4082: SNESComputeFunction_Matlab - Calls the function that has been set with
4083: SNESSetFunctionMatlab().
4085: Collective on SNES
4087: Input Parameters:
4088: + snes - the SNES context
4089: - x - input vector
4091: Output Parameter:
4092: . y - function vector, as set by SNESSetFunction()
4094: Notes:
4095: SNESComputeFunction() is typically used within nonlinear solvers
4096: implementations, so most users would not generally call this routine
4097: themselves.
4099: Level: developer
4101: .keywords: SNES, nonlinear, compute, function
4103: .seealso: SNESSetFunction(), SNESGetFunction()
4104: */
4105: PetscErrorCode SNESComputeFunction_Matlab(SNES snes,Vec x,Vec y, void *ctx)
4106: {
4107: PetscErrorCode ierr;
4108: SNESMatlabContext *sctx = (SNESMatlabContext *)ctx;
4109: int nlhs = 1,nrhs = 5;
4110: mxArray *plhs[1],*prhs[5];
4111: long long int lx = 0,ly = 0,ls = 0;
4120: /* call Matlab function in ctx with arguments x and y */
4122: PetscMemcpy(&ls,&snes,sizeof(snes));
4123: PetscMemcpy(&lx,&x,sizeof(x));
4124: PetscMemcpy(&ly,&y,sizeof(x));
4125: prhs[0] = mxCreateDoubleScalar((double)ls);
4126: prhs[1] = mxCreateDoubleScalar((double)lx);
4127: prhs[2] = mxCreateDoubleScalar((double)ly);
4128: prhs[3] = mxCreateString(sctx->funcname);
4129: prhs[4] = sctx->ctx;
4130: mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeFunctionInternal");
4131: mxGetScalar(plhs[0]);
4132: mxDestroyArray(prhs[0]);
4133: mxDestroyArray(prhs[1]);
4134: mxDestroyArray(prhs[2]);
4135: mxDestroyArray(prhs[3]);
4136: mxDestroyArray(plhs[0]);
4137: return(0);
4138: }
4143: /*
4144: SNESSetFunctionMatlab - Sets the function evaluation routine and function
4145: vector for use by the SNES routines in solving systems of nonlinear
4146: equations from MATLAB. Here the function is a string containing the name of a MATLAB function
4148: Logically Collective on SNES
4150: Input Parameters:
4151: + snes - the SNES context
4152: . r - vector to store function value
4153: - func - function evaluation routine
4155: Calling sequence of func:
4156: $ func (SNES snes,Vec x,Vec f,void *ctx);
4159: Notes:
4160: The Newton-like methods typically solve linear systems of the form
4161: $ f'(x) x = -f(x),
4162: where f'(x) denotes the Jacobian matrix and f(x) is the function.
4164: Level: beginner
4166: .keywords: SNES, nonlinear, set, function
4168: .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction()
4169: */
4170: PetscErrorCode SNESSetFunctionMatlab(SNES snes,Vec r,const char *func,mxArray *ctx)
4171: {
4172: PetscErrorCode ierr;
4173: SNESMatlabContext *sctx;
4176: /* currently sctx is memory bleed */
4177: PetscMalloc(sizeof(SNESMatlabContext),&sctx);
4178: PetscStrallocpy(func,&sctx->funcname);
4179: /*
4180: This should work, but it doesn't
4181: sctx->ctx = ctx;
4182: mexMakeArrayPersistent(sctx->ctx);
4183: */
4184: sctx->ctx = mxDuplicateArray(ctx);
4185: SNESSetFunction(snes,r,SNESComputeFunction_Matlab,sctx);
4186: return(0);
4187: }
4191: /*
4192: SNESComputeJacobian_Matlab - Calls the function that has been set with
4193: SNESSetJacobianMatlab().
4195: Collective on SNES
4197: Input Parameters:
4198: + snes - the SNES context
4199: . x - input vector
4200: . A, B - the matrices
4201: - ctx - user context
4203: Output Parameter:
4204: . flag - structure of the matrix
4206: Level: developer
4208: .keywords: SNES, nonlinear, compute, function
4210: .seealso: SNESSetFunction(), SNESGetFunction()
4211: @*/
4212: PetscErrorCode SNESComputeJacobian_Matlab(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *flag, void *ctx)
4213: {
4214: PetscErrorCode ierr;
4215: SNESMatlabContext *sctx = (SNESMatlabContext *)ctx;
4216: int nlhs = 2,nrhs = 6;
4217: mxArray *plhs[2],*prhs[6];
4218: long long int lx = 0,lA = 0,ls = 0, lB = 0;
4219:
4224: /* call Matlab function in ctx with arguments x and y */
4226: PetscMemcpy(&ls,&snes,sizeof(snes));
4227: PetscMemcpy(&lx,&x,sizeof(x));
4228: PetscMemcpy(&lA,A,sizeof(x));
4229: PetscMemcpy(&lB,B,sizeof(x));
4230: prhs[0] = mxCreateDoubleScalar((double)ls);
4231: prhs[1] = mxCreateDoubleScalar((double)lx);
4232: prhs[2] = mxCreateDoubleScalar((double)lA);
4233: prhs[3] = mxCreateDoubleScalar((double)lB);
4234: prhs[4] = mxCreateString(sctx->funcname);
4235: prhs[5] = sctx->ctx;
4236: mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeJacobianInternal");
4237: mxGetScalar(plhs[0]);
4238: *flag = (MatStructure) mxGetScalar(plhs[1]);
4239: mxDestroyArray(prhs[0]);
4240: mxDestroyArray(prhs[1]);
4241: mxDestroyArray(prhs[2]);
4242: mxDestroyArray(prhs[3]);
4243: mxDestroyArray(prhs[4]);
4244: mxDestroyArray(plhs[0]);
4245: mxDestroyArray(plhs[1]);
4246: return(0);
4247: }
4252: /*
4253: SNESSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices
4254: vector for use by the SNES routines in solving systems of nonlinear
4255: equations from MATLAB. Here the function is a string containing the name of a MATLAB function
4257: Logically Collective on SNES
4259: Input Parameters:
4260: + snes - the SNES context
4261: . A,B - Jacobian matrices
4262: . func - function evaluation routine
4263: - ctx - user context
4265: Calling sequence of func:
4266: $ flag = func (SNES snes,Vec x,Mat A,Mat B,void *ctx);
4269: Level: developer
4271: .keywords: SNES, nonlinear, set, function
4273: .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction()
4274: */
4275: PetscErrorCode SNESSetJacobianMatlab(SNES snes,Mat A,Mat B,const char *func,mxArray *ctx)
4276: {
4277: PetscErrorCode ierr;
4278: SNESMatlabContext *sctx;
4281: /* currently sctx is memory bleed */
4282: PetscMalloc(sizeof(SNESMatlabContext),&sctx);
4283: PetscStrallocpy(func,&sctx->funcname);
4284: /*
4285: This should work, but it doesn't
4286: sctx->ctx = ctx;
4287: mexMakeArrayPersistent(sctx->ctx);
4288: */
4289: sctx->ctx = mxDuplicateArray(ctx);
4290: SNESSetJacobian(snes,A,B,SNESComputeJacobian_Matlab,sctx);
4291: return(0);
4292: }
4296: /*
4297: SNESMonitor_Matlab - Calls the function that has been set with SNESMonitorSetMatlab().
4299: Collective on SNES
4301: .seealso: SNESSetFunction(), SNESGetFunction()
4302: @*/
4303: PetscErrorCode SNESMonitor_Matlab(SNES snes,PetscInt it, PetscReal fnorm, void *ctx)
4304: {
4305: PetscErrorCode ierr;
4306: SNESMatlabContext *sctx = (SNESMatlabContext *)ctx;
4307: int nlhs = 1,nrhs = 6;
4308: mxArray *plhs[1],*prhs[6];
4309: long long int lx = 0,ls = 0;
4310: Vec x=snes->vec_sol;
4311:
4315: PetscMemcpy(&ls,&snes,sizeof(snes));
4316: PetscMemcpy(&lx,&x,sizeof(x));
4317: prhs[0] = mxCreateDoubleScalar((double)ls);
4318: prhs[1] = mxCreateDoubleScalar((double)it);
4319: prhs[2] = mxCreateDoubleScalar((double)fnorm);
4320: prhs[3] = mxCreateDoubleScalar((double)lx);
4321: prhs[4] = mxCreateString(sctx->funcname);
4322: prhs[5] = sctx->ctx;
4323: mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESMonitorInternal");
4324: mxGetScalar(plhs[0]);
4325: mxDestroyArray(prhs[0]);
4326: mxDestroyArray(prhs[1]);
4327: mxDestroyArray(prhs[2]);
4328: mxDestroyArray(prhs[3]);
4329: mxDestroyArray(prhs[4]);
4330: mxDestroyArray(plhs[0]);
4331: return(0);
4332: }
4337: /*
4338: SNESMonitorSetMatlab - Sets the monitor function from MATLAB
4340: Level: developer
4342: .keywords: SNES, nonlinear, set, function
4344: .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction()
4345: */
4346: PetscErrorCode SNESMonitorSetMatlab(SNES snes,const char *func,mxArray *ctx)
4347: {
4348: PetscErrorCode ierr;
4349: SNESMatlabContext *sctx;
4352: /* currently sctx is memory bleed */
4353: PetscMalloc(sizeof(SNESMatlabContext),&sctx);
4354: PetscStrallocpy(func,&sctx->funcname);
4355: /*
4356: This should work, but it doesn't
4357: sctx->ctx = ctx;
4358: mexMakeArrayPersistent(sctx->ctx);
4359: */
4360: sctx->ctx = mxDuplicateArray(ctx);
4361: SNESMonitorSet(snes,SNESMonitor_Matlab,sctx,PETSC_NULL);
4362: return(0);
4363: }
4365: #endif