Actual source code: snes.c
petsc-3.3-p7 2013-05-11
2: #include <petsc-private/snesimpl.h> /*I "petscsnes.h" I*/
3: #include <petscdmshell.h> /*I "petscdmshell.h" I*/
4: #include <petscsys.h> /*I "petscsys.h" I*/
6: PetscBool SNESRegisterAllCalled = PETSC_FALSE;
7: PetscFList SNESList = PETSC_NULL;
9: /* Logging support */
10: PetscClassId SNES_CLASSID;
11: PetscLogEvent SNES_Solve, SNES_FunctionEval, SNES_JacobianEval, SNES_GSEval;
15: /*@
16: SNESSetErrorIfNotConverged - Causes SNESSolve() to generate an error if the solver has not converged.
18: Logically Collective on SNES
20: Input Parameters:
21: + snes - iterative context obtained from SNESCreate()
22: - flg - PETSC_TRUE indicates you want the error generated
24: Options database keys:
25: . -snes_error_if_not_converged : this takes an optional truth value (0/1/no/yes/true/false)
27: Level: intermediate
29: Notes:
30: Normally PETSc continues if a linear solver fails to converge, you can call SNESGetConvergedReason() after a SNESSolve()
31: to determine if it has converged.
33: .keywords: SNES, set, initial guess, nonzero
35: .seealso: SNESGetErrorIfNotConverged(), KSPGetErrorIfNotConverged(), KSPSetErrorIFNotConverged()
36: @*/
37: PetscErrorCode SNESSetErrorIfNotConverged(SNES snes,PetscBool flg)
38: {
42: snes->errorifnotconverged = flg;
44: return(0);
45: }
49: /*@
50: SNESGetErrorIfNotConverged - Will SNESSolve() generate an error if the solver does not converge?
52: Not Collective
54: Input Parameter:
55: . snes - iterative context obtained from SNESCreate()
57: Output Parameter:
58: . flag - PETSC_TRUE if it will generate an error, else PETSC_FALSE
60: Level: intermediate
62: .keywords: SNES, set, initial guess, nonzero
64: .seealso: SNESSetErrorIfNotConverged(), KSPGetErrorIfNotConverged(), KSPSetErrorIFNotConverged()
65: @*/
66: PetscErrorCode SNESGetErrorIfNotConverged(SNES snes,PetscBool *flag)
67: {
71: *flag = snes->errorifnotconverged;
72: return(0);
73: }
77: /*@
78: SNESSetFunctionDomainError - tells SNES that the input vector to your FormFunction is not
79: in the functions domain. For example, negative pressure.
81: Logically Collective on SNES
83: Input Parameters:
84: . snes - the SNES context
86: Level: advanced
88: .keywords: SNES, view
90: .seealso: SNESCreate(), SNESSetFunction()
91: @*/
92: PetscErrorCode SNESSetFunctionDomainError(SNES snes)
93: {
96: snes->domainerror = PETSC_TRUE;
97: return(0);
98: }
103: /*@
104: SNESGetFunctionDomainError - Gets the status of the domain error after a call to SNESComputeFunction;
106: Logically Collective on SNES
108: Input Parameters:
109: . snes - the SNES context
111: Output Parameters:
112: . domainerror Set to PETSC_TRUE if there's a domain error; PETSC_FALSE otherwise.
114: Level: advanced
116: .keywords: SNES, view
118: .seealso: SNESSetFunctionDomainError, SNESComputeFunction()
119: @*/
120: PetscErrorCode SNESGetFunctionDomainError(SNES snes, PetscBool *domainerror)
121: {
125: *domainerror = snes->domainerror;
126: return(0);
127: }
132: /*@C
133: SNESView - Prints the SNES data structure.
135: Collective on SNES
137: Input Parameters:
138: + SNES - the SNES context
139: - viewer - visualization context
141: Options Database Key:
142: . -snes_view - Calls SNESView() at end of SNESSolve()
144: Notes:
145: The available visualization contexts include
146: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
147: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
148: output where only the first processor opens
149: the file. All other processors send their
150: data to the first processor to print.
152: The user can open an alternative visualization context with
153: PetscViewerASCIIOpen() - output to a specified file.
155: Level: beginner
157: .keywords: SNES, view
159: .seealso: PetscViewerASCIIOpen()
160: @*/
161: PetscErrorCode SNESView(SNES snes,PetscViewer viewer)
162: {
163: SNESKSPEW *kctx;
164: PetscErrorCode ierr;
165: KSP ksp;
166: SNESLineSearch linesearch;
167: PetscBool iascii,isstring;
171: if (!viewer) {
172: PetscViewerASCIIGetStdout(((PetscObject)snes)->comm,&viewer);
173: }
177: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
178: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
179: if (iascii) {
180: PetscObjectPrintClassNamePrefixType((PetscObject)snes,viewer,"SNES Object");
181: if (snes->ops->view) {
182: PetscViewerASCIIPushTab(viewer);
183: (*snes->ops->view)(snes,viewer);
184: PetscViewerASCIIPopTab(viewer);
185: }
186: PetscViewerASCIIPrintf(viewer," maximum iterations=%D, maximum function evaluations=%D\n",snes->max_its,snes->max_funcs);
187: PetscViewerASCIIPrintf(viewer," tolerances: relative=%G, absolute=%G, solution=%G\n",
188: snes->rtol,snes->abstol,snes->stol);
189: PetscViewerASCIIPrintf(viewer," total number of linear solver iterations=%D\n",snes->linear_its);
190: PetscViewerASCIIPrintf(viewer," total number of function evaluations=%D\n",snes->nfuncs);
191: if (snes->gridsequence) {
192: PetscViewerASCIIPrintf(viewer," total number of grid sequence refinements=%D\n",snes->gridsequence);
193: }
194: if (snes->ksp_ewconv) {
195: kctx = (SNESKSPEW *)snes->kspconvctx;
196: if (kctx) {
197: PetscViewerASCIIPrintf(viewer," Eisenstat-Walker computation of KSP relative tolerance (version %D)\n",kctx->version);
198: PetscViewerASCIIPrintf(viewer," rtol_0=%G, rtol_max=%G, threshold=%G\n",kctx->rtol_0,kctx->rtol_max,kctx->threshold);
199: PetscViewerASCIIPrintf(viewer," gamma=%G, alpha=%G, alpha2=%G\n",kctx->gamma,kctx->alpha,kctx->alpha2);
200: }
201: }
202: if (snes->lagpreconditioner == -1) {
203: PetscViewerASCIIPrintf(viewer," Preconditioned is never rebuilt\n");
204: } else if (snes->lagpreconditioner > 1) {
205: PetscViewerASCIIPrintf(viewer," Preconditioned is rebuilt every %D new Jacobians\n",snes->lagpreconditioner);
206: }
207: if (snes->lagjacobian == -1) {
208: PetscViewerASCIIPrintf(viewer," Jacobian is never rebuilt\n");
209: } else if (snes->lagjacobian > 1) {
210: PetscViewerASCIIPrintf(viewer," Jacobian is rebuilt every %D SNES iterations\n",snes->lagjacobian);
211: }
212: } else if (isstring) {
213: const char *type;
214: SNESGetType(snes,&type);
215: PetscViewerStringSPrintf(viewer," %-3.3s",type);
216: }
217: if (snes->pc && snes->usespc) {
218: PetscViewerASCIIPushTab(viewer);
219: SNESView(snes->pc, viewer);
220: PetscViewerASCIIPopTab(viewer);
221: }
222: if (snes->usesksp) {
223: SNESGetKSP(snes,&ksp);
224: PetscViewerASCIIPushTab(viewer);
225: KSPView(ksp,viewer);
226: PetscViewerASCIIPopTab(viewer);
227: }
228: if (snes->linesearch) {
229: PetscViewerASCIIPushTab(viewer);
230: SNESGetSNESLineSearch(snes, &linesearch);
231: SNESLineSearchView(linesearch, viewer);
232: PetscViewerASCIIPopTab(viewer);
233: }
234: return(0);
235: }
237: /*
238: We retain a list of functions that also take SNES command
239: line options. These are called at the end SNESSetFromOptions()
240: */
241: #define MAXSETFROMOPTIONS 5
242: static PetscInt numberofsetfromoptions;
243: static PetscErrorCode (*othersetfromoptions[MAXSETFROMOPTIONS])(SNES);
247: /*@C
248: SNESAddOptionsChecker - Adds an additional function to check for SNES options.
250: Not Collective
252: Input Parameter:
253: . snescheck - function that checks for options
255: Level: developer
257: .seealso: SNESSetFromOptions()
258: @*/
259: PetscErrorCode SNESAddOptionsChecker(PetscErrorCode (*snescheck)(SNES))
260: {
262: if (numberofsetfromoptions >= MAXSETFROMOPTIONS) {
263: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Too many options checkers, only %D allowed", MAXSETFROMOPTIONS);
264: }
265: othersetfromoptions[numberofsetfromoptions++] = snescheck;
266: return(0);
267: }
269: extern PetscErrorCode SNESDefaultMatrixFreeCreate2(SNES,Vec,Mat*);
273: static PetscErrorCode SNESSetUpMatrixFree_Private(SNES snes, PetscBool hasOperator, PetscInt version)
274: {
275: Mat J;
276: KSP ksp;
277: PC pc;
278: PetscBool match;
284: if(!snes->vec_func && (snes->jacobian || snes->jacobian_pre)) {
285: Mat A = snes->jacobian, B = snes->jacobian_pre;
286: MatGetVecs(A ? A : B, PETSC_NULL,&snes->vec_func);
287: }
289: if (version == 1) {
290: MatCreateSNESMF(snes,&J);
291: MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);
292: MatSetFromOptions(J);
293: } else if (version == 2) {
294: if (!snes->vec_func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"SNESSetFunction() must be called first");
295: #if !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128)
296: SNESDefaultMatrixFreeCreate2(snes,snes->vec_func,&J);
297: #else
298: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "matrix-free operator rutines (version 2)");
299: #endif
300: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "matrix-free operator rutines, only version 1 and 2");
301:
302: PetscInfo1(snes,"Setting default matrix-free operator routines (version %D)\n", version);
303: if (hasOperator) {
304: /* This version replaces the user provided Jacobian matrix with a
305: matrix-free version but still employs the user-provided preconditioner matrix. */
306: SNESSetJacobian(snes,J,0,0,0);
307: } else {
308: /* This version replaces both the user-provided Jacobian and the user-
309: provided preconditioner matrix with the default matrix free version. */
310: void *functx;
311: SNESGetFunction(snes,PETSC_NULL,PETSC_NULL,&functx);
312: SNESSetJacobian(snes,J,J,MatMFFDComputeJacobian,functx);
313: /* Force no preconditioner */
314: SNESGetKSP(snes,&ksp);
315: KSPGetPC(ksp,&pc);
316: PetscObjectTypeCompare((PetscObject)pc,PCSHELL,&match);
317: if (!match) {
318: PetscInfo(snes,"Setting default matrix-free preconditioner routines\nThat is no preconditioner is being used\n");
319: PCSetType(pc,PCNONE);
320: }
321: }
322: MatDestroy(&J);
323: return(0);
324: }
328: static PetscErrorCode DMRestrictHook_SNESVecSol(DM dmfine,Mat Restrict,Vec Rscale,Mat Inject,DM dmcoarse,void *ctx)
329: {
330: SNES snes = (SNES)ctx;
332: Vec Xfine,Xfine_named = PETSC_NULL,Xcoarse;
335: if (PetscLogPrintInfo) {
336: PetscInt finelevel,coarselevel,fineclevel,coarseclevel;
337: DMGetRefineLevel(dmfine,&finelevel);
338: DMGetCoarsenLevel(dmfine,&fineclevel);
339: DMGetRefineLevel(dmcoarse,&coarselevel);
340: DMGetCoarsenLevel(dmcoarse,&coarseclevel);
341: PetscInfo4(dmfine,"Restricting SNES solution vector from level %D-%D to level %D-%D\n",finelevel,fineclevel,coarselevel,coarseclevel);
342: }
343: if (dmfine == snes->dm) Xfine = snes->vec_sol;
344: else {
345: DMGetNamedGlobalVector(dmfine,"SNESVecSol",&Xfine_named);
346: Xfine = Xfine_named;
347: }
348: DMGetNamedGlobalVector(dmcoarse,"SNESVecSol",&Xcoarse);
349: MatRestrict(Restrict,Xfine,Xcoarse);
350: VecPointwiseMult(Xcoarse,Xcoarse,Rscale);
351: DMRestoreNamedGlobalVector(dmcoarse,"SNESVecSol",&Xcoarse);
352: if (Xfine_named) {DMRestoreNamedGlobalVector(dmfine,"SNESVecSol",&Xfine_named);}
353: return(0);
354: }
358: static PetscErrorCode DMCoarsenHook_SNESVecSol(DM dm,DM dmc,void *ctx)
359: {
363: DMCoarsenHookAdd(dmc,DMCoarsenHook_SNESVecSol,DMRestrictHook_SNESVecSol,ctx);
364: return(0);
365: }
369: /* This may be called to rediscretize the operator on levels of linear multigrid. The DM shuffle is so the user can
370: * safely call SNESGetDM() in their residual evaluation routine. */
371: static PetscErrorCode KSPComputeOperators_SNES(KSP ksp,Mat A,Mat B,MatStructure *mstruct,void *ctx)
372: {
373: SNES snes = (SNES)ctx;
375: Mat Asave = A,Bsave = B;
376: Vec X,Xnamed = PETSC_NULL;
377: DM dmsave;
380: dmsave = snes->dm;
381: KSPGetDM(ksp,&snes->dm);
382: if (dmsave == snes->dm) X = snes->vec_sol; /* We are on the finest level */
383: else { /* We are on a coarser level, this vec was initialized using a DM restrict hook */
384: DMGetNamedGlobalVector(snes->dm,"SNESVecSol",&Xnamed);
385: X = Xnamed;
386: }
387: SNESComputeJacobian(snes,X,&A,&B,mstruct);
388: if (A != Asave || B != Bsave) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_SUP,"No support for changing matrices at this time");
389: if (Xnamed) {
390: DMRestoreNamedGlobalVector(snes->dm,"SNESVecSol",&Xnamed);
391: }
392: snes->dm = dmsave;
393: return(0);
394: }
398: /*@
399: SNESSetUpMatrices - ensures that matrices are available for SNES, to be called by SNESSetUp_XXX()
401: Collective
403: Input Arguments:
404: . snes - snes to configure
406: Level: developer
408: .seealso: SNESSetUp()
409: @*/
410: PetscErrorCode SNESSetUpMatrices(SNES snes)
411: {
413: DM dm;
414: SNESDM sdm;
417: SNESGetDM(snes,&dm);
418: DMSNESGetContext(dm,&sdm);
419: if (!sdm->computejacobian) {
420: SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_PLIB,"SNESDM not improperly configured");
421: } else if (!snes->jacobian && sdm->computejacobian == MatMFFDComputeJacobian) {
422: Mat J;
423: void *functx;
424: MatCreateSNESMF(snes,&J);
425: MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);
426: MatSetFromOptions(J);
427: SNESGetFunction(snes,PETSC_NULL,PETSC_NULL,&functx);
428: SNESSetJacobian(snes,J,J,MatMFFDComputeJacobian,functx);
429: MatDestroy(&J);
430: } else if (snes->mf_operator && !snes->jacobian_pre && !snes->jacobian) {
431: Mat J,B;
432: MatCreateSNESMF(snes,&J);
433: MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);
434: MatSetFromOptions(J);
435: DMCreateMatrix(snes->dm,MATAIJ,&B);
436: /* sdm->computejacobian was already set to reach here */
437: SNESSetJacobian(snes,J,B,PETSC_NULL,PETSC_NULL);
438: MatDestroy(&J);
439: MatDestroy(&B);
440: } else if (!snes->jacobian_pre) {
441: Mat J,B;
442: J = snes->jacobian;
443: DMCreateMatrix(snes->dm,MATAIJ,&B);
444: SNESSetJacobian(snes,J?J:B,B,PETSC_NULL,PETSC_NULL);
445: MatDestroy(&B);
446: }
447: {
448: KSP ksp;
449: SNESGetKSP(snes,&ksp);
450: KSPSetComputeOperators(ksp,KSPComputeOperators_SNES,snes);
451: DMCoarsenHookAdd(snes->dm,DMCoarsenHook_SNESVecSol,DMRestrictHook_SNESVecSol,snes);
452: }
453: return(0);
454: }
458: /*@
459: SNESSetFromOptions - Sets various SNES and KSP parameters from user options.
461: Collective on SNES
463: Input Parameter:
464: . snes - the SNES context
466: Options Database Keys:
467: + -snes_type <type> - ls, tr, ngmres, ncg, richardson, qn, vi, fas
468: . -snes_stol - convergence tolerance in terms of the norm
469: of the change in the solution between steps
470: . -snes_atol <abstol> - absolute tolerance of residual norm
471: . -snes_rtol <rtol> - relative decrease in tolerance norm from initial
472: . -snes_max_it <max_it> - maximum number of iterations
473: . -snes_max_funcs <max_funcs> - maximum number of function evaluations
474: . -snes_max_fail <max_fail> - maximum number of line search failures allowed before stopping, default is none
475: . -snes_max_linear_solve_fail - number of linear solver failures before SNESSolve() stops
476: . -snes_lag_preconditioner <lag> - how often preconditioner is rebuilt (use -1 to never rebuild)
477: . -snes_lag_jacobian <lag> - how often Jacobian is rebuilt (use -1 to never rebuild)
478: . -snes_trtol <trtol> - trust region tolerance
479: . -snes_no_convergence_test - skip convergence test in nonlinear
480: solver; hence iterations will continue until max_it
481: or some other criterion is reached. Saves expense
482: of convergence test
483: . -snes_monitor <optional filename> - prints residual norm at each iteration. if no
484: filename given prints to stdout
485: . -snes_monitor_solution - plots solution at each iteration
486: . -snes_monitor_residual - plots residual (not its norm) at each iteration
487: . -snes_monitor_solution_update - plots update to solution at each iteration
488: . -snes_monitor_draw - plots residual norm at each iteration
489: . -snes_fd - use finite differences to compute Jacobian; very slow, only for testing
490: . -snes_mf_ksp_monitor - if using matrix-free multiply then print h at each KSP iteration
491: - -snes_converged_reason - print the reason for convergence/divergence after each solve
493: Options Database for Eisenstat-Walker method:
494: + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence
495: . -snes_ksp_ew_version ver - version of Eisenstat-Walker method
496: . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0
497: . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax
498: . -snes_ksp_ew_gamma <gamma> - Sets gamma
499: . -snes_ksp_ew_alpha <alpha> - Sets alpha
500: . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2
501: - -snes_ksp_ew_threshold <threshold> - Sets threshold
503: Notes:
504: To see all options, run your program with the -help option or consult
505: the <A href="../../docs/manual.pdf#nameddest=ch_snes">SNES chapter of the users manual</A>.
507: Level: beginner
509: .keywords: SNES, nonlinear, set, options, database
511: .seealso: SNESSetOptionsPrefix()
512: @*/
513: PetscErrorCode SNESSetFromOptions(SNES snes)
514: {
515: PetscBool flg,mf,mf_operator,pcset;
516: PetscInt i,indx,lag,mf_version,grids;
517: MatStructure matflag;
518: const char *deft = SNESLS;
519: const char *convtests[] = {"default","skip"};
520: SNESKSPEW *kctx = NULL;
521: char type[256], monfilename[PETSC_MAX_PATH_LEN];
522: PetscViewer monviewer;
523: PetscErrorCode ierr;
524: const char *optionsprefix;
529: if (!SNESRegisterAllCalled) {SNESRegisterAll(PETSC_NULL);}
530: PetscObjectOptionsBegin((PetscObject)snes);
531: if (((PetscObject)snes)->type_name) { deft = ((PetscObject)snes)->type_name; }
532: PetscOptionsList("-snes_type","Nonlinear solver method","SNESSetType",SNESList,deft,type,256,&flg);
533: if (flg) {
534: SNESSetType(snes,type);
535: } else if (!((PetscObject)snes)->type_name) {
536: SNESSetType(snes,deft);
537: }
538: /* not used here, but called so will go into help messaage */
539: PetscOptionsName("-snes_view","Print detailed information on solver used","SNESView",0);
541: PetscOptionsReal("-snes_stol","Stop if step length less than","SNESSetTolerances",snes->stol,&snes->stol,0);
542: PetscOptionsReal("-snes_atol","Stop if function norm less than","SNESSetTolerances",snes->abstol,&snes->abstol,0);
544: PetscOptionsReal("-snes_rtol","Stop if decrease in function norm less than","SNESSetTolerances",snes->rtol,&snes->rtol,0);
545: PetscOptionsInt("-snes_max_it","Maximum iterations","SNESSetTolerances",snes->max_its,&snes->max_its,PETSC_NULL);
546: PetscOptionsInt("-snes_max_funcs","Maximum function evaluations","SNESSetTolerances",snes->max_funcs,&snes->max_funcs,PETSC_NULL);
547: PetscOptionsInt("-snes_max_fail","Maximum nonlinear step failures","SNESSetMaxNonlinearStepFailures",snes->maxFailures,&snes->maxFailures,PETSC_NULL);
548: PetscOptionsInt("-snes_max_linear_solve_fail","Maximum failures in linear solves allowed","SNESSetMaxLinearSolveFailures",snes->maxLinearSolveFailures,&snes->maxLinearSolveFailures,PETSC_NULL);
549: PetscOptionsBool("-snes_error_if_not_converged","Generate error if solver does not converge","SNESSetErrorIfNotConverged",snes->errorifnotconverged,&snes->errorifnotconverged,PETSC_NULL);
551: PetscOptionsInt("-snes_lag_preconditioner","How often to rebuild preconditioner","SNESSetLagPreconditioner",snes->lagpreconditioner,&lag,&flg);
552: if (flg) {
553: SNESSetLagPreconditioner(snes,lag);
554: }
555: PetscOptionsInt("-snes_lag_jacobian","How often to rebuild Jacobian","SNESSetLagJacobian",snes->lagjacobian,&lag,&flg);
556: if (flg) {
557: SNESSetLagJacobian(snes,lag);
558: }
559: PetscOptionsInt("-snes_grid_sequence","Use grid sequencing to generate initial guess","SNESSetGridSequence",snes->gridsequence,&grids,&flg);
560: if (flg) {
561: SNESSetGridSequence(snes,grids);
562: }
564: PetscOptionsEList("-snes_convergence_test","Convergence test","SNESSetConvergenceTest",convtests,2,"default",&indx,&flg);
565: if (flg) {
566: switch (indx) {
567: case 0: SNESSetConvergenceTest(snes,SNESDefaultConverged,PETSC_NULL,PETSC_NULL); break;
568: case 1: SNESSetConvergenceTest(snes,SNESSkipConverged,PETSC_NULL,PETSC_NULL); break;
569: }
570: }
572: PetscOptionsBool("-snes_converged_reason","Print reason for converged or diverged","SNESSolve",snes->printreason,&snes->printreason,PETSC_NULL);
574: PetscOptionsEList("-snes_norm_type","SNES Norm type","SNESSetNormType",SNESNormTypes,5,"function",&indx,&flg);
575: if (flg) { SNESSetNormType(snes,(SNESNormType)indx); }
577: kctx = (SNESKSPEW *)snes->kspconvctx;
579: PetscOptionsBool("-snes_ksp_ew","Use Eisentat-Walker linear system convergence test","SNESKSPSetUseEW",snes->ksp_ewconv,&snes->ksp_ewconv,PETSC_NULL);
581: PetscOptionsInt("-snes_ksp_ew_version","Version 1, 2 or 3","SNESKSPSetParametersEW",kctx->version,&kctx->version,0);
582: PetscOptionsReal("-snes_ksp_ew_rtol0","0 <= rtol0 < 1","SNESKSPSetParametersEW",kctx->rtol_0,&kctx->rtol_0,0);
583: PetscOptionsReal("-snes_ksp_ew_rtolmax","0 <= rtolmax < 1","SNESKSPSetParametersEW",kctx->rtol_max,&kctx->rtol_max,0);
584: PetscOptionsReal("-snes_ksp_ew_gamma","0 <= gamma <= 1","SNESKSPSetParametersEW",kctx->gamma,&kctx->gamma,0);
585: PetscOptionsReal("-snes_ksp_ew_alpha","1 < alpha <= 2","SNESKSPSetParametersEW",kctx->alpha,&kctx->alpha,0);
586: PetscOptionsReal("-snes_ksp_ew_alpha2","alpha2","SNESKSPSetParametersEW",kctx->alpha2,&kctx->alpha2,0);
587: PetscOptionsReal("-snes_ksp_ew_threshold","0 < threshold < 1","SNESKSPSetParametersEW",kctx->threshold,&kctx->threshold,0);
589: flg = PETSC_FALSE;
590: PetscOptionsBool("-snes_monitor_cancel","Remove all monitors","SNESMonitorCancel",flg,&flg,PETSC_NULL);
591: if (flg) {SNESMonitorCancel(snes);}
593: PetscOptionsString("-snes_monitor","Monitor norm of function","SNESMonitorSet","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);
594: if (flg) {
595: PetscViewerASCIIOpen(((PetscObject)snes)->comm,monfilename,&monviewer);
596: SNESMonitorSet(snes,SNESMonitorDefault,monviewer,(PetscErrorCode (*)(void**))PetscViewerDestroy);
597: }
599: PetscOptionsString("-snes_monitor_range","Monitor range of elements of function","SNESMonitorSet","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);
600: if (flg) {
601: SNESMonitorSet(snes,SNESMonitorRange,0,0);
602: }
604: PetscOptionsString("-snes_ratiomonitor","Monitor ratios of norms of function","SNESMonitorSetRatio","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);
605: if (flg) {
606: PetscViewerASCIIOpen(((PetscObject)snes)->comm,monfilename,&monviewer);
607: SNESMonitorSetRatio(snes,monviewer);
608: }
610: PetscOptionsString("-snes_monitor_short","Monitor norm of function (fewer digits)","SNESMonitorSet","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);
611: if (flg) {
612: PetscViewerASCIIOpen(((PetscObject)snes)->comm,monfilename,&monviewer);
613: SNESMonitorSet(snes,SNESMonitorDefaultShort,monviewer,(PetscErrorCode (*)(void**))PetscViewerDestroy);
614: }
616: PetscOptionsString("-snes_monitor_python","Use Python function","SNESMonitorSet",0,monfilename,PETSC_MAX_PATH_LEN,&flg);
617: if (flg) {PetscPythonMonitorSet((PetscObject)snes,monfilename);}
619: flg = PETSC_FALSE;
620: PetscOptionsBool("-snes_monitor_solution","Plot solution at each iteration","SNESMonitorSolution",flg,&flg,PETSC_NULL);
621: if (flg) {SNESMonitorSet(snes,SNESMonitorSolution,0,0);}
622: flg = PETSC_FALSE;
623: PetscOptionsBool("-snes_monitor_solution_update","Plot correction at each iteration","SNESMonitorSolutionUpdate",flg,&flg,PETSC_NULL);
624: if (flg) {SNESMonitorSet(snes,SNESMonitorSolutionUpdate,0,0);}
625: flg = PETSC_FALSE;
626: PetscOptionsBool("-snes_monitor_residual","Plot residual at each iteration","SNESMonitorResidual",flg,&flg,PETSC_NULL);
627: if (flg) {SNESMonitorSet(snes,SNESMonitorResidual,0,0);}
628: flg = PETSC_FALSE;
629: PetscOptionsBool("-snes_monitor_draw","Plot function norm at each iteration","SNESMonitorLG",flg,&flg,PETSC_NULL);
630: if (flg) {SNESMonitorSet(snes,SNESMonitorLG,PETSC_NULL,PETSC_NULL);}
631: flg = PETSC_FALSE;
632: PetscOptionsBool("-snes_monitor_range_draw","Plot function range at each iteration","SNESMonitorLG",flg,&flg,PETSC_NULL);
633: if (flg) {SNESMonitorSet(snes,SNESMonitorLGRange,PETSC_NULL,PETSC_NULL);}
635: flg = PETSC_FALSE;
636: PetscOptionsBool("-snes_fd","Use finite differences (slow) to compute Jacobian","SNESDefaultComputeJacobian",flg,&flg,PETSC_NULL);
637: if (flg) {
638: void *functx;
639: SNESGetFunction(snes,PETSC_NULL,PETSC_NULL,&functx);
640: SNESSetJacobian(snes,snes->jacobian,snes->jacobian_pre,SNESDefaultComputeJacobian,functx);
641: PetscInfo(snes,"Setting default finite difference Jacobian matrix\n");
642: }
644: mf = mf_operator = PETSC_FALSE;
645: flg = PETSC_FALSE;
646: PetscOptionsBool("-snes_mf_operator","Use a Matrix-Free Jacobian with user-provided preconditioner matrix","MatCreateSNESMF",PETSC_FALSE,&mf_operator,&flg);
647: if (flg && mf_operator) {
648: snes->mf_operator = PETSC_TRUE;
649: mf = PETSC_TRUE;
650: }
651: flg = PETSC_FALSE;
652: PetscOptionsBool("-snes_mf","Use a Matrix-Free Jacobian with no preconditioner matrix","MatCreateSNESMF",PETSC_FALSE,&mf,&flg);
653: if (!flg && mf_operator) mf = PETSC_TRUE;
654: mf_version = 1;
655: PetscOptionsInt("-snes_mf_version","Matrix-Free routines version 1 or 2","None",mf_version,&mf_version,0);
658: /* GS Options */
659: PetscOptionsInt("-snes_gs_sweeps","Number of sweeps of GS to apply","SNESComputeGS",snes->gssweeps,&snes->gssweeps,PETSC_NULL);
661: for(i = 0; i < numberofsetfromoptions; i++) {
662: (*othersetfromoptions[i])(snes);
663: }
665: if (snes->ops->setfromoptions) {
666: (*snes->ops->setfromoptions)(snes);
667: }
669: /* process any options handlers added with PetscObjectAddOptionsHandler() */
670: PetscObjectProcessOptionsHandlers((PetscObject)snes);
671: PetscOptionsEnd();
673: if (mf) { SNESSetUpMatrixFree_Private(snes, mf_operator, mf_version); }
675: if (!snes->ksp) {SNESGetKSP(snes,&snes->ksp);}
676: KSPGetOperators(snes->ksp,PETSC_NULL,PETSC_NULL,&matflag);
677: KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,matflag);
678: KSPSetFromOptions(snes->ksp);
680: if (!snes->linesearch) {
681: SNESGetSNESLineSearch(snes, &snes->linesearch);
682: }
683: SNESLineSearchSetFromOptions(snes->linesearch);
685: /* if someone has set the SNES PC type, create it. */
686: SNESGetOptionsPrefix(snes, &optionsprefix);
687: PetscOptionsHasName(optionsprefix, "-npc_snes_type", &pcset);
688: if (pcset && (!snes->pc)) {
689: SNESGetPC(snes, &snes->pc);
690: }
691: return(0);
692: }
696: /*@
697: SNESSetComputeApplicationContext - Sets an optional function to compute a user-defined context for
698: the nonlinear solvers.
700: Logically Collective on SNES
702: Input Parameters:
703: + snes - the SNES context
704: . compute - function to compute the context
705: - destroy - function to destroy the context
707: Level: intermediate
709: .keywords: SNES, nonlinear, set, application, context
711: .seealso: SNESGetApplicationContext(), SNESSetComputeApplicationContext(), SNESGetApplicationContext()
712: @*/
713: PetscErrorCode SNESSetComputeApplicationContext(SNES snes,PetscErrorCode (*compute)(SNES,void**),PetscErrorCode (*destroy)(void**))
714: {
717: snes->ops->usercompute = compute;
718: snes->ops->userdestroy = destroy;
719: return(0);
720: }
724: /*@
725: SNESSetApplicationContext - Sets the optional user-defined context for
726: the nonlinear solvers.
728: Logically Collective on SNES
730: Input Parameters:
731: + snes - the SNES context
732: - usrP - optional user context
734: Level: intermediate
736: .keywords: SNES, nonlinear, set, application, context
738: .seealso: SNESGetApplicationContext()
739: @*/
740: PetscErrorCode SNESSetApplicationContext(SNES snes,void *usrP)
741: {
743: KSP ksp;
747: SNESGetKSP(snes,&ksp);
748: KSPSetApplicationContext(ksp,usrP);
749: snes->user = usrP;
750: return(0);
751: }
755: /*@
756: SNESGetApplicationContext - Gets the user-defined context for the
757: nonlinear solvers.
759: Not Collective
761: Input Parameter:
762: . snes - SNES context
764: Output Parameter:
765: . usrP - user context
767: Level: intermediate
769: .keywords: SNES, nonlinear, get, application, context
771: .seealso: SNESSetApplicationContext()
772: @*/
773: PetscErrorCode SNESGetApplicationContext(SNES snes,void *usrP)
774: {
777: *(void**)usrP = snes->user;
778: return(0);
779: }
783: /*@
784: SNESGetIterationNumber - Gets the number of nonlinear iterations completed
785: at this time.
787: Not Collective
789: Input Parameter:
790: . snes - SNES context
792: Output Parameter:
793: . iter - iteration number
795: Notes:
796: For example, during the computation of iteration 2 this would return 1.
798: This is useful for using lagged Jacobians (where one does not recompute the
799: Jacobian at each SNES iteration). For example, the code
800: .vb
801: SNESGetIterationNumber(snes,&it);
802: if (!(it % 2)) {
803: [compute Jacobian here]
804: }
805: .ve
806: can be used in your ComputeJacobian() function to cause the Jacobian to be
807: recomputed every second SNES iteration.
809: Level: intermediate
811: .keywords: SNES, nonlinear, get, iteration, number,
813: .seealso: SNESGetFunctionNorm(), SNESGetLinearSolveIterations()
814: @*/
815: PetscErrorCode SNESGetIterationNumber(SNES snes,PetscInt* iter)
816: {
820: *iter = snes->iter;
821: return(0);
822: }
826: /*@
827: SNESSetIterationNumber - Sets the current iteration number.
829: Not Collective
831: Input Parameter:
832: . snes - SNES context
833: . iter - iteration number
835: Level: developer
837: .keywords: SNES, nonlinear, set, iteration, number,
839: .seealso: SNESGetFunctionNorm(), SNESGetLinearSolveIterations()
840: @*/
841: PetscErrorCode SNESSetIterationNumber(SNES snes,PetscInt iter)
842: {
847: PetscObjectTakeAccess(snes);
848: snes->iter = iter;
849: PetscObjectGrantAccess(snes);
850: return(0);
851: }
855: /*@
856: SNESGetFunctionNorm - Gets the norm of the current function that was set
857: with SNESSSetFunction().
859: Collective on SNES
861: Input Parameter:
862: . snes - SNES context
864: Output Parameter:
865: . fnorm - 2-norm of function
867: Level: intermediate
869: .keywords: SNES, nonlinear, get, function, norm
871: .seealso: SNESGetFunction(), SNESGetIterationNumber(), SNESGetLinearSolveIterations()
872: @*/
873: PetscErrorCode SNESGetFunctionNorm(SNES snes,PetscReal *fnorm)
874: {
878: *fnorm = snes->norm;
879: return(0);
880: }
885: /*@
886: SNESSetFunctionNorm - Sets the 2-norm of the current function computed using VecNorm().
888: Collective on SNES
890: Input Parameter:
891: . snes - SNES context
892: . fnorm - 2-norm of function
894: Level: developer
896: .keywords: SNES, nonlinear, set, function, norm
898: .seealso: SNESSetFunction(), SNESSetIterationNumber(), VecNorm().
899: @*/
900: PetscErrorCode SNESSetFunctionNorm(SNES snes,PetscReal fnorm)
901: {
907: PetscObjectTakeAccess(snes);
908: snes->norm = fnorm;
909: PetscObjectGrantAccess(snes);
910: return(0);
911: }
915: /*@
916: SNESGetNonlinearStepFailures - Gets the number of unsuccessful steps
917: attempted by the nonlinear solver.
919: Not Collective
921: Input Parameter:
922: . snes - SNES context
924: Output Parameter:
925: . nfails - number of unsuccessful steps attempted
927: Notes:
928: This counter is reset to zero for each successive call to SNESSolve().
930: Level: intermediate
932: .keywords: SNES, nonlinear, get, number, unsuccessful, steps
934: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures(),
935: SNESSetMaxNonlinearStepFailures(), SNESGetMaxNonlinearStepFailures()
936: @*/
937: PetscErrorCode SNESGetNonlinearStepFailures(SNES snes,PetscInt* nfails)
938: {
942: *nfails = snes->numFailures;
943: return(0);
944: }
948: /*@
949: SNESSetMaxNonlinearStepFailures - Sets the maximum number of unsuccessful steps
950: attempted by the nonlinear solver before it gives up.
952: Not Collective
954: Input Parameters:
955: + snes - SNES context
956: - maxFails - maximum of unsuccessful steps
958: Level: intermediate
960: .keywords: SNES, nonlinear, set, maximum, unsuccessful, steps
962: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures(),
963: SNESGetMaxNonlinearStepFailures(), SNESGetNonlinearStepFailures()
964: @*/
965: PetscErrorCode SNESSetMaxNonlinearStepFailures(SNES snes, PetscInt maxFails)
966: {
969: snes->maxFailures = maxFails;
970: return(0);
971: }
975: /*@
976: SNESGetMaxNonlinearStepFailures - Gets the maximum number of unsuccessful steps
977: attempted by the nonlinear solver before it gives up.
979: Not Collective
981: Input Parameter:
982: . snes - SNES context
984: Output Parameter:
985: . maxFails - maximum of unsuccessful steps
987: Level: intermediate
989: .keywords: SNES, nonlinear, get, maximum, unsuccessful, steps
991: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures(),
992: SNESSetMaxNonlinearStepFailures(), SNESGetNonlinearStepFailures()
993:
994: @*/
995: PetscErrorCode SNESGetMaxNonlinearStepFailures(SNES snes, PetscInt *maxFails)
996: {
1000: *maxFails = snes->maxFailures;
1001: return(0);
1002: }
1006: /*@
1007: SNESGetNumberFunctionEvals - Gets the number of user provided function evaluations
1008: done by SNES.
1010: Not Collective
1012: Input Parameter:
1013: . snes - SNES context
1015: Output Parameter:
1016: . nfuncs - number of evaluations
1018: Level: intermediate
1020: .keywords: SNES, nonlinear, get, maximum, unsuccessful, steps
1022: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures()
1023: @*/
1024: PetscErrorCode SNESGetNumberFunctionEvals(SNES snes, PetscInt *nfuncs)
1025: {
1029: *nfuncs = snes->nfuncs;
1030: return(0);
1031: }
1035: /*@
1036: SNESGetLinearSolveFailures - Gets the number of failed (non-converged)
1037: linear solvers.
1039: Not Collective
1041: Input Parameter:
1042: . snes - SNES context
1044: Output Parameter:
1045: . nfails - number of failed solves
1047: Notes:
1048: This counter is reset to zero for each successive call to SNESSolve().
1050: Level: intermediate
1052: .keywords: SNES, nonlinear, get, number, unsuccessful, steps
1054: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures()
1055: @*/
1056: PetscErrorCode SNESGetLinearSolveFailures(SNES snes,PetscInt* nfails)
1057: {
1061: *nfails = snes->numLinearSolveFailures;
1062: return(0);
1063: }
1067: /*@
1068: SNESSetMaxLinearSolveFailures - the number of failed linear solve attempts
1069: allowed before SNES returns with a diverged reason of SNES_DIVERGED_LINEAR_SOLVE
1071: Logically Collective on SNES
1073: Input Parameters:
1074: + snes - SNES context
1075: - maxFails - maximum allowed linear solve failures
1077: Level: intermediate
1079: Notes: By default this is 0; that is SNES returns on the first failed linear solve
1081: .keywords: SNES, nonlinear, set, maximum, unsuccessful, steps
1083: .seealso: SNESGetLinearSolveFailures(), SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations()
1084: @*/
1085: PetscErrorCode SNESSetMaxLinearSolveFailures(SNES snes, PetscInt maxFails)
1086: {
1090: snes->maxLinearSolveFailures = maxFails;
1091: return(0);
1092: }
1096: /*@
1097: SNESGetMaxLinearSolveFailures - gets the maximum number of linear solve failures that
1098: are allowed before SNES terminates
1100: Not Collective
1102: Input Parameter:
1103: . snes - SNES context
1105: Output Parameter:
1106: . maxFails - maximum of unsuccessful solves allowed
1108: Level: intermediate
1110: Notes: By default this is 1; that is SNES returns on the first failed linear solve
1112: .keywords: SNES, nonlinear, get, maximum, unsuccessful, steps
1114: .seealso: SNESGetLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(),
1115: @*/
1116: PetscErrorCode SNESGetMaxLinearSolveFailures(SNES snes, PetscInt *maxFails)
1117: {
1121: *maxFails = snes->maxLinearSolveFailures;
1122: return(0);
1123: }
1127: /*@
1128: SNESGetLinearSolveIterations - Gets the total number of linear iterations
1129: used by the nonlinear solver.
1131: Not Collective
1133: Input Parameter:
1134: . snes - SNES context
1136: Output Parameter:
1137: . lits - number of linear iterations
1139: Notes:
1140: This counter is reset to zero for each successive call to SNESSolve().
1142: Level: intermediate
1144: .keywords: SNES, nonlinear, get, number, linear, iterations
1146: .seealso: SNESGetIterationNumber(), SNESGetFunctionNorm(), SNESGetLinearSolveFailures(), SNESGetMaxLinearSolveFailures()
1147: @*/
1148: PetscErrorCode SNESGetLinearSolveIterations(SNES snes,PetscInt* lits)
1149: {
1153: *lits = snes->linear_its;
1154: return(0);
1155: }
1159: /*@
1160: SNESGetKSP - Returns the KSP context for a SNES solver.
1162: Not Collective, but if SNES object is parallel, then KSP object is parallel
1164: Input Parameter:
1165: . snes - the SNES context
1167: Output Parameter:
1168: . ksp - the KSP context
1170: Notes:
1171: The user can then directly manipulate the KSP context to set various
1172: options, etc. Likewise, the user can then extract and manipulate the
1173: PC contexts as well.
1175: Level: beginner
1177: .keywords: SNES, nonlinear, get, KSP, context
1179: .seealso: KSPGetPC(), SNESCreate(), KSPCreate(), SNESSetKSP()
1180: @*/
1181: PetscErrorCode SNESGetKSP(SNES snes,KSP *ksp)
1182: {
1189: if (!snes->ksp) {
1190: KSPCreate(((PetscObject)snes)->comm,&snes->ksp);
1191: PetscObjectIncrementTabLevel((PetscObject)snes->ksp,(PetscObject)snes,1);
1192: PetscLogObjectParent(snes,snes->ksp);
1193: }
1194: *ksp = snes->ksp;
1195: return(0);
1196: }
1200: /*@
1201: SNESSetKSP - Sets a KSP context for the SNES object to use
1203: Not Collective, but the SNES and KSP objects must live on the same MPI_Comm
1205: Input Parameters:
1206: + snes - the SNES context
1207: - ksp - the KSP context
1209: Notes:
1210: The SNES object already has its KSP object, you can obtain with SNESGetKSP()
1211: so this routine is rarely needed.
1213: The KSP object that is already in the SNES object has its reference count
1214: decreased by one.
1216: Level: developer
1218: .keywords: SNES, nonlinear, get, KSP, context
1220: .seealso: KSPGetPC(), SNESCreate(), KSPCreate(), SNESSetKSP()
1221: @*/
1222: PetscErrorCode SNESSetKSP(SNES snes,KSP ksp)
1223: {
1230: PetscObjectReference((PetscObject)ksp);
1231: if (snes->ksp) {PetscObjectDereference((PetscObject)snes->ksp);}
1232: snes->ksp = ksp;
1233: return(0);
1234: }
1236: #if 0
1239: static PetscErrorCode SNESPublish_Petsc(PetscObject obj)
1240: {
1242: return(0);
1243: }
1244: #endif
1246: /* -----------------------------------------------------------*/
1249: /*@
1250: SNESCreate - Creates a nonlinear solver context.
1252: Collective on MPI_Comm
1254: Input Parameters:
1255: . comm - MPI communicator
1257: Output Parameter:
1258: . outsnes - the new SNES context
1260: Options Database Keys:
1261: + -snes_mf - Activates default matrix-free Jacobian-vector products,
1262: and no preconditioning matrix
1263: . -snes_mf_operator - Activates default matrix-free Jacobian-vector
1264: products, and a user-provided preconditioning matrix
1265: as set by SNESSetJacobian()
1266: - -snes_fd - Uses (slow!) finite differences to compute Jacobian
1268: Level: beginner
1270: .keywords: SNES, nonlinear, create, context
1272: .seealso: SNESSolve(), SNESDestroy(), SNES, SNESSetLagPreconditioner()
1274: @*/
1275: PetscErrorCode SNESCreate(MPI_Comm comm,SNES *outsnes)
1276: {
1277: PetscErrorCode ierr;
1278: SNES snes;
1279: SNESKSPEW *kctx;
1283: *outsnes = PETSC_NULL;
1284: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
1285: SNESInitializePackage(PETSC_NULL);
1286: #endif
1288: PetscHeaderCreate(snes,_p_SNES,struct _SNESOps,SNES_CLASSID,0,"SNES","Nonlinear solver","SNES",comm,SNESDestroy,SNESView);
1290: snes->ops->converged = SNESDefaultConverged;
1291: snes->usesksp = PETSC_TRUE;
1292: snes->tolerancesset = PETSC_FALSE;
1293: snes->max_its = 50;
1294: snes->max_funcs = 10000;
1295: snes->norm = 0.0;
1296: snes->normtype = SNES_NORM_FUNCTION;
1297: snes->rtol = 1.e-8;
1298: snes->ttol = 0.0;
1299: snes->abstol = 1.e-50;
1300: snes->stol = 1.e-8;
1301: snes->deltatol = 1.e-12;
1302: snes->nfuncs = 0;
1303: snes->numFailures = 0;
1304: snes->maxFailures = 1;
1305: snes->linear_its = 0;
1306: snes->lagjacobian = 1;
1307: snes->lagpreconditioner = 1;
1308: snes->numbermonitors = 0;
1309: snes->data = 0;
1310: snes->setupcalled = PETSC_FALSE;
1311: snes->ksp_ewconv = PETSC_FALSE;
1312: snes->nwork = 0;
1313: snes->work = 0;
1314: snes->nvwork = 0;
1315: snes->vwork = 0;
1316: snes->conv_hist_len = 0;
1317: snes->conv_hist_max = 0;
1318: snes->conv_hist = PETSC_NULL;
1319: snes->conv_hist_its = PETSC_NULL;
1320: snes->conv_hist_reset = PETSC_TRUE;
1321: snes->vec_func_init_set = PETSC_FALSE;
1322: snes->norm_init = 0.;
1323: snes->norm_init_set = PETSC_FALSE;
1324: snes->reason = SNES_CONVERGED_ITERATING;
1325: snes->gssweeps = 1;
1327: snes->numLinearSolveFailures = 0;
1328: snes->maxLinearSolveFailures = 1;
1330: /* Create context to compute Eisenstat-Walker relative tolerance for KSP */
1331: PetscNewLog(snes,SNESKSPEW,&kctx);
1332: snes->kspconvctx = (void*)kctx;
1333: kctx->version = 2;
1334: kctx->rtol_0 = .3; /* Eisenstat and Walker suggest rtol_0=.5, but
1335: this was too large for some test cases */
1336: kctx->rtol_last = 0.0;
1337: kctx->rtol_max = .9;
1338: kctx->gamma = 1.0;
1339: kctx->alpha = .5*(1.0 + PetscSqrtReal(5.0));
1340: kctx->alpha2 = kctx->alpha;
1341: kctx->threshold = .1;
1342: kctx->lresid_last = 0.0;
1343: kctx->norm_last = 0.0;
1345: *outsnes = snes;
1346: return(0);
1347: }
1351: /*@C
1352: SNESSetFunction - Sets the function evaluation routine and function
1353: vector for use by the SNES routines in solving systems of nonlinear
1354: equations.
1356: Logically Collective on SNES
1358: Input Parameters:
1359: + snes - the SNES context
1360: . r - vector to store function value
1361: . func - function evaluation routine
1362: - ctx - [optional] user-defined context for private data for the
1363: function evaluation routine (may be PETSC_NULL)
1365: Calling sequence of func:
1366: $ func (SNES snes,Vec x,Vec f,void *ctx);
1368: + snes - the SNES context
1369: . x - state at which to evaluate residual
1370: . f - vector to put residual
1371: - ctx - optional user-defined function context
1373: Notes:
1374: The Newton-like methods typically solve linear systems of the form
1375: $ f'(x) x = -f(x),
1376: where f'(x) denotes the Jacobian matrix and f(x) is the function.
1378: Level: beginner
1380: .keywords: SNES, nonlinear, set, function
1382: .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetPicard()
1383: @*/
1384: PetscErrorCode SNESSetFunction(SNES snes,Vec r,PetscErrorCode (*func)(SNES,Vec,Vec,void*),void *ctx)
1385: {
1387: DM dm;
1391: if (r) {
1394: PetscObjectReference((PetscObject)r);
1395: VecDestroy(&snes->vec_func);
1396: snes->vec_func = r;
1397: }
1398: SNESGetDM(snes,&dm);
1399: DMSNESSetFunction(dm,func,ctx);
1400: return(0);
1401: }
1406: /*@C
1407: SNESSetInitialFunction - Sets the function vector to be used as the
1408: function norm at the initialization of the method. In some
1409: instances, the user has precomputed the function before calling
1410: SNESSolve. This function allows one to avoid a redundant call
1411: to SNESComputeFunction in that case.
1413: Logically Collective on SNES
1415: Input Parameters:
1416: + snes - the SNES context
1417: - f - vector to store function value
1419: Notes:
1420: This should not be modified during the solution procedure.
1422: This is used extensively in the SNESFAS hierarchy and in nonlinear preconditioning.
1424: Level: developer
1426: .keywords: SNES, nonlinear, set, function
1428: .seealso: SNESSetFunction(), SNESComputeFunction(), SNESSetInitialFunctionNorm()
1429: @*/
1430: PetscErrorCode SNESSetInitialFunction(SNES snes, Vec f)
1431: {
1433: Vec vec_func;
1439: SNESGetFunction(snes,&vec_func,PETSC_NULL,PETSC_NULL);
1440: VecCopy(f, vec_func);
1441: snes->vec_func_init_set = PETSC_TRUE;
1442: return(0);
1443: }
1448: /*@C
1449: SNESSetInitialFunctionNorm - Sets the function norm to be used as the function norm
1450: at the initialization of the method. In some instances, the user has precomputed
1451: the function and its norm before calling SNESSolve. This function allows one to
1452: avoid a redundant call to SNESComputeFunction() and VecNorm() in that case.
1454: Logically Collective on SNES
1456: Input Parameters:
1457: + snes - the SNES context
1458: - fnorm - the norm of F as set by SNESSetInitialFunction()
1460: This is used extensively in the SNESFAS hierarchy and in nonlinear preconditioning.
1462: Level: developer
1464: .keywords: SNES, nonlinear, set, function, norm
1466: .seealso: SNESSetFunction(), SNESComputeFunction(), SNESSetInitialFunction()
1467: @*/
1468: PetscErrorCode SNESSetInitialFunctionNorm(SNES snes, PetscReal fnorm)
1469: {
1472: snes->norm_init = fnorm;
1473: snes->norm_init_set = PETSC_TRUE;
1474: return(0);
1475: }
1479: /*@
1480: SNESSetNormType - Sets the SNESNormType used in covergence and monitoring
1481: of the SNES method.
1483: Logically Collective on SNES
1485: Input Parameters:
1486: + snes - the SNES context
1487: - normtype - the type of the norm used
1489: Notes:
1490: Only certain SNES methods support certain SNESNormTypes. Most require evaluation
1491: of the nonlinear function and the taking of its norm at every iteration to
1492: even ensure convergence at all. However, methods such as custom Gauss-Seidel methods
1493: (SNESGS) and the like do not require the norm of the function to be computed, and therfore
1494: may either be monitored for convergence or not. As these are often used as nonlinear
1495: preconditioners, monitoring the norm of their error is not a useful enterprise within
1496: their solution.
1498: Level: developer
1500: .keywords: SNES, nonlinear, set, function, norm, type
1502: .seealso: SNESGetNormType(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormType
1503: @*/
1504: PetscErrorCode SNESSetNormType(SNES snes, SNESNormType normtype)
1505: {
1508: snes->normtype = normtype;
1509: return(0);
1510: }
1515: /*@
1516: SNESGetNormType - Gets the SNESNormType used in covergence and monitoring
1517: of the SNES method.
1519: Logically Collective on SNES
1521: Input Parameters:
1522: + snes - the SNES context
1523: - normtype - the type of the norm used
1525: Level: advanced
1527: .keywords: SNES, nonlinear, set, function, norm, type
1529: .seealso: SNESSetNormType(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormType
1530: @*/
1531: PetscErrorCode SNESGetNormType(SNES snes, SNESNormType *normtype)
1532: {
1535: *normtype = snes->normtype;
1536: return(0);
1537: }
1541: /*@C
1542: SNESSetGS - Sets the user nonlinear Gauss-Seidel routine for
1543: use with composed nonlinear solvers.
1545: Input Parameters:
1546: + snes - the SNES context
1547: . gsfunc - function evaluation routine
1548: - ctx - [optional] user-defined context for private data for the
1549: smoother evaluation routine (may be PETSC_NULL)
1551: Calling sequence of func:
1552: $ func (SNES snes,Vec x,Vec b,void *ctx);
1554: + X - solution vector
1555: . B - RHS vector
1556: - ctx - optional user-defined Gauss-Seidel context
1558: Notes:
1559: The GS routines are used by the composed nonlinear solver to generate
1560: a problem appropriate update to the solution, particularly FAS.
1562: Level: intermediate
1564: .keywords: SNES, nonlinear, set, Gauss-Seidel
1566: .seealso: SNESGetFunction(), SNESComputeGS()
1567: @*/
1568: PetscErrorCode SNESSetGS(SNES snes,PetscErrorCode (*gsfunc)(SNES,Vec,Vec,void*),void *ctx)
1569: {
1571: DM dm;
1575: SNESGetDM(snes,&dm);
1576: DMSNESSetGS(dm,gsfunc,ctx);
1577: return(0);
1578: }
1582: /*@
1583: SNESSetGSSweeps - Sets the number of sweeps of GS to use.
1585: Input Parameters:
1586: + snes - the SNES context
1587: - sweeps - the number of sweeps of GS to perform.
1589: Level: intermediate
1591: .keywords: SNES, nonlinear, set, Gauss-Siedel
1593: .seealso: SNESSetGS(), SNESGetGS(), SNESSetPC(), SNESGetGSSweeps()
1594: @*/
1596: PetscErrorCode SNESSetGSSweeps(SNES snes, PetscInt sweeps) {
1598: snes->gssweeps = sweeps;
1599: return(0);
1600: }
1605: /*@
1606: SNESGetGSSweeps - Gets the number of sweeps GS will use.
1608: Input Parameters:
1609: . snes - the SNES context
1611: Output Parameters:
1612: . sweeps - the number of sweeps of GS to perform.
1614: Level: intermediate
1616: .keywords: SNES, nonlinear, set, Gauss-Siedel
1618: .seealso: SNESSetGS(), SNESGetGS(), SNESSetPC(), SNESSetGSSweeps()
1619: @*/
1620: PetscErrorCode SNESGetGSSweeps(SNES snes, PetscInt * sweeps) {
1622: *sweeps = snes->gssweeps;
1623: return(0);
1624: }
1628: PetscErrorCode SNESPicardComputeFunction(SNES snes,Vec x,Vec f,void *ctx)
1629: {
1631: DM dm;
1632: SNESDM sdm;
1635: SNESGetDM(snes,&dm);
1636: DMSNESGetContext(dm,&sdm);
1637: /* A(x)*x - b(x) */
1638: if (sdm->computepfunction) {
1639: (*sdm->computepfunction)(snes,x,f,sdm->pctx);
1640: } else if (snes->dm) {
1641: DMComputeFunction(snes->dm,x,f);
1642: } else {
1643: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetPicard() or SNESSetDM() before SNESPicardComputeFunction to provide Picard function.");
1644: }
1646: if (sdm->computepjacobian) {
1647: (*sdm->computepjacobian)(snes,x,&snes->jacobian,&snes->jacobian_pre,&snes->matstruct,sdm->pctx);
1648: } else if (snes->dm) {
1649: DMComputeJacobian(snes->dm,x,snes->jacobian,snes->jacobian_pre,&snes->matstruct);
1650: } else {
1651: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetPicard() or SNESSetDM() before SNESPicardComputeFunction to provide Picard matrix.");
1652: }
1654: VecView(x,PETSC_VIEWER_BINARY_WORLD);
1655: VecView(f,PETSC_VIEWER_BINARY_WORLD);
1656: VecScale(f,-1.0);
1657: MatMultAdd(snes->jacobian_pre,x,f,f);
1658: return(0);
1659: }
1663: PetscErrorCode SNESPicardComputeJacobian(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx)
1664: {
1666: /* the jacobian matrix should be pre-filled in SNESPicardComputeFunction */
1667: *flag = snes->matstruct;
1668: return(0);
1669: }
1673: /*@C
1674: SNESSetPicard - Use SNES to solve the semilinear-system A(x) x = b(x) via a Picard type iteration (Picard linearization)
1676: Logically Collective on SNES
1678: Input Parameters:
1679: + snes - the SNES context
1680: . r - vector to store function value
1681: . func - function evaluation routine
1682: . 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)
1683: . mat - matrix to store A
1684: . mfunc - function to compute matrix value
1685: - ctx - [optional] user-defined context for private data for the
1686: function evaluation routine (may be PETSC_NULL)
1688: Calling sequence of func:
1689: $ func (SNES snes,Vec x,Vec f,void *ctx);
1691: + f - function vector
1692: - ctx - optional user-defined function context
1694: Calling sequence of mfunc:
1695: $ mfunc (SNES snes,Vec x,Mat *jmat,Mat *mat,int *flag,void *ctx);
1697: + x - input vector
1698: . 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(),
1699: normally just pass mat in this location
1700: . mat - form A(x) matrix
1701: . flag - flag indicating information about the preconditioner matrix
1702: structure (same as flag in KSPSetOperators()), one of SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER
1703: - ctx - [optional] user-defined Jacobian context
1705: Notes:
1706: One can call SNESSetPicard() or SNESSetFunction() (and possibly SNESSetJacobian()) but cannot call both
1708: $ 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}
1709: $ Note that when an exact solver is used this corresponds to the "classic" Picard A(x^{n}) x^{n+1} = b(x^{n}) iteration.
1711: Run with -snes_mf_operator to solve the system with Newton's method using A(x^{n}) to construct the preconditioner.
1713: We implement the defect correction form of the Picard iteration because it converges much more generally when inexact linear solvers are used then
1714: the direct Picard iteration A(x^n) x^{n+1} = b(x^n)
1716: 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
1717: 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
1718: different please contact us at petsc-dev@mcs.anl.gov and we'll have an entirely new argument :-).
1720: Level: beginner
1722: .keywords: SNES, nonlinear, set, function
1724: .seealso: SNESGetFunction(), SNESSetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESGetPicard(), SNESLineSearchPreCheckPicard()
1725: @*/
1726: 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)
1727: {
1729: DM dm;
1733: SNESGetDM(snes, &dm);
1734: DMSNESSetPicard(dm,func,mfunc,ctx);
1735: SNESSetFunction(snes,r,SNESPicardComputeFunction,ctx);
1736: SNESSetJacobian(snes,jmat,mat,SNESPicardComputeJacobian,ctx);
1737: return(0);
1738: }
1743: /*@C
1744: SNESGetPicard - Returns the context for the Picard iteration
1746: Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet.
1748: Input Parameter:
1749: . snes - the SNES context
1751: Output Parameter:
1752: + r - the function (or PETSC_NULL)
1753: . func - the function (or PETSC_NULL)
1754: . jmat - the picard matrix (or PETSC_NULL)
1755: . mat - the picard preconditioner matrix (or PETSC_NULL)
1756: . mfunc - the function for matrix evaluation (or PETSC_NULL)
1757: - ctx - the function context (or PETSC_NULL)
1759: Level: advanced
1761: .keywords: SNES, nonlinear, get, function
1763: .seealso: SNESSetPicard, SNESGetFunction, SNESGetJacobian, SNESGetDM
1764: @*/
1765: PetscErrorCode SNESGetPicard(SNES snes,Vec *r,PetscErrorCode (**func)(SNES,Vec,Vec,void*),Mat *jmat, Mat *mat, PetscErrorCode (**mfunc)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void **ctx)
1766: {
1768: DM dm;
1772: SNESGetFunction(snes,r,PETSC_NULL,PETSC_NULL);
1773: SNESGetJacobian(snes,jmat,mat,PETSC_NULL,PETSC_NULL);
1774: SNESGetDM(snes,&dm);
1775: DMSNESGetPicard(dm,func,mfunc,ctx);
1776: return(0);
1777: }
1781: /*@C
1782: SNESSetComputeInitialGuess - Sets a routine used to compute an initial guess for the problem
1784: Logically Collective on SNES
1786: Input Parameters:
1787: + snes - the SNES context
1788: . func - function evaluation routine
1789: - ctx - [optional] user-defined context for private data for the
1790: function evaluation routine (may be PETSC_NULL)
1792: Calling sequence of func:
1793: $ func (SNES snes,Vec x,void *ctx);
1795: . f - function vector
1796: - ctx - optional user-defined function context
1798: Level: intermediate
1800: .keywords: SNES, nonlinear, set, function
1802: .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian()
1803: @*/
1804: PetscErrorCode SNESSetComputeInitialGuess(SNES snes,PetscErrorCode (*func)(SNES,Vec,void*),void *ctx)
1805: {
1808: if (func) snes->ops->computeinitialguess = func;
1809: if (ctx) snes->initialguessP = ctx;
1810: return(0);
1811: }
1813: /* --------------------------------------------------------------- */
1816: /*@C
1817: SNESGetRhs - Gets the vector for solving F(x) = rhs. If rhs is not set
1818: it assumes a zero right hand side.
1820: Logically Collective on SNES
1822: Input Parameter:
1823: . snes - the SNES context
1825: Output Parameter:
1826: . rhs - the right hand side vector or PETSC_NULL if the right hand side vector is null
1828: Level: intermediate
1830: .keywords: SNES, nonlinear, get, function, right hand side
1832: .seealso: SNESGetSolution(), SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction()
1833: @*/
1834: PetscErrorCode SNESGetRhs(SNES snes,Vec *rhs)
1835: {
1839: *rhs = snes->vec_rhs;
1840: return(0);
1841: }
1845: /*@
1846: SNESComputeFunction - Calls the function that has been set with
1847: SNESSetFunction().
1849: Collective on SNES
1851: Input Parameters:
1852: + snes - the SNES context
1853: - x - input vector
1855: Output Parameter:
1856: . y - function vector, as set by SNESSetFunction()
1858: Notes:
1859: SNESComputeFunction() is typically used within nonlinear solvers
1860: implementations, so most users would not generally call this routine
1861: themselves.
1863: Level: developer
1865: .keywords: SNES, nonlinear, compute, function
1867: .seealso: SNESSetFunction(), SNESGetFunction()
1868: @*/
1869: PetscErrorCode SNESComputeFunction(SNES snes,Vec x,Vec y)
1870: {
1872: DM dm;
1873: SNESDM sdm;
1881: VecValidValues(x,2,PETSC_TRUE);
1883: SNESGetDM(snes,&dm);
1884: DMSNESGetContext(dm,&sdm);
1885: PetscLogEventBegin(SNES_FunctionEval,snes,x,y,0);
1886: if (sdm->computefunction) {
1887: PetscStackPush("SNES user function");
1888: (*sdm->computefunction)(snes,x,y,sdm->functionctx);
1889: PetscStackPop;
1890: } else if (snes->dm) {
1891: DMComputeFunction(snes->dm,x,y);
1892: } else if (snes->vec_rhs) {
1893: MatMult(snes->jacobian, x, y);
1894: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetFunction() or SNESSetDM() before SNESComputeFunction(), likely called from SNESSolve().");
1895: if (snes->vec_rhs) {
1896: VecAXPY(y,-1.0,snes->vec_rhs);
1897: }
1898: snes->nfuncs++;
1899: PetscLogEventEnd(SNES_FunctionEval,snes,x,y,0);
1900: VecValidValues(y,3,PETSC_FALSE);
1901: return(0);
1902: }
1906: /*@
1907: SNESComputeGS - Calls the Gauss-Seidel function that has been set with
1908: SNESSetGS().
1910: Collective on SNES
1912: Input Parameters:
1913: + snes - the SNES context
1914: . x - input vector
1915: - b - rhs vector
1917: Output Parameter:
1918: . x - new solution vector
1920: Notes:
1921: SNESComputeGS() is typically used within composed nonlinear solver
1922: implementations, so most users would not generally call this routine
1923: themselves.
1925: Level: developer
1927: .keywords: SNES, nonlinear, compute, function
1929: .seealso: SNESSetGS(), SNESComputeFunction()
1930: @*/
1931: PetscErrorCode SNESComputeGS(SNES snes,Vec b,Vec x)
1932: {
1934: PetscInt i;
1935: DM dm;
1936: SNESDM sdm;
1944: if (b) VecValidValues(b,2,PETSC_TRUE);
1945: PetscLogEventBegin(SNES_GSEval,snes,x,b,0);
1946: SNESGetDM(snes,&dm);
1947: DMSNESGetContext(dm,&sdm);
1948: if (sdm->computegs) {
1949: for (i = 0; i < snes->gssweeps; i++) {
1950: PetscStackPush("SNES user GS");
1951: (*sdm->computegs)(snes,x,b,sdm->gsctx);
1952: PetscStackPop;
1953: }
1954: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetGS() before SNESComputeGS(), likely called from SNESSolve().");
1955: PetscLogEventEnd(SNES_GSEval,snes,x,b,0);
1956: VecValidValues(x,3,PETSC_FALSE);
1957: return(0);
1958: }
1963: /*@
1964: SNESComputeJacobian - Computes the Jacobian matrix that has been
1965: set with SNESSetJacobian().
1967: Collective on SNES and Mat
1969: Input Parameters:
1970: + snes - the SNES context
1971: - x - input vector
1973: Output Parameters:
1974: + A - Jacobian matrix
1975: . B - optional preconditioning matrix
1976: - flag - flag indicating matrix structure (one of, SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER)
1978: Options Database Keys:
1979: + -snes_lag_preconditioner <lag>
1980: . -snes_lag_jacobian <lag>
1981: . -snes_compare_explicit - Compare the computed Jacobian to the finite difference Jacobian and output the differences
1982: . -snes_compare_explicit_draw - Compare the computed Jacobian to the finite difference Jacobian and draw the result
1983: . -snes_compare_explicit_contour - Compare the computed Jacobian to the finite difference Jacobian and draw a contour plot with the result
1984: . -snes_compare_operator - Make the comparison options above use the operator instead of the preconditioning matrix
1985: . -snes_compare_coloring - Compute the finite differece Jacobian using coloring and display norms of difference
1986: . -snes_compare_coloring_display - Compute the finite differece Jacobian using coloring and display verbose differences
1987: . -snes_compare_coloring_threshold - Display only those matrix entries that differ by more than a given threshold
1988: . -snes_compare_coloring_threshold_atol - Absolute tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold
1989: . -snes_compare_coloring_threshold_rtol - Relative tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold
1990: . -snes_compare_coloring_draw - Compute the finite differece Jacobian using coloring and draw differences
1991: - -snes_compare_coloring_draw_contour - Compute the finite differece Jacobian using coloring and show contours of matrices and differences
1994: Notes:
1995: Most users should not need to explicitly call this routine, as it
1996: is used internally within the nonlinear solvers.
1998: See KSPSetOperators() for important information about setting the
1999: flag parameter.
2001: Level: developer
2003: .keywords: SNES, compute, Jacobian, matrix
2005: .seealso: SNESSetJacobian(), KSPSetOperators(), MatStructure, SNESSetLagPreconditioner(), SNESSetLagJacobian()
2006: @*/
2007: PetscErrorCode SNESComputeJacobian(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *flg)
2008: {
2010: PetscBool flag;
2011: DM dm;
2012: SNESDM sdm;
2019: VecValidValues(X,2,PETSC_TRUE);
2020: SNESGetDM(snes,&dm);
2021: DMSNESGetContext(dm,&sdm);
2022: if (!sdm->computejacobian) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_USER,"Must call SNESSetJacobian(), DMSNESSetJacobian(), DMDASNESSetJacobianLocal(), etc");
2024: /* make sure that MatAssemblyBegin/End() is called on A matrix if it is matrix free */
2026: if (snes->lagjacobian == -2) {
2027: snes->lagjacobian = -1;
2028: PetscInfo(snes,"Recomputing Jacobian/preconditioner because lag is -2 (means compute Jacobian, but then never again) \n");
2029: } else if (snes->lagjacobian == -1) {
2030: *flg = SAME_PRECONDITIONER;
2031: PetscInfo(snes,"Reusing Jacobian/preconditioner because lag is -1\n");
2032: PetscObjectTypeCompare((PetscObject)*A,MATMFFD,&flag);
2033: if (flag) {
2034: MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
2035: MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
2036: }
2037: return(0);
2038: } else if (snes->lagjacobian > 1 && snes->iter % snes->lagjacobian) {
2039: *flg = SAME_PRECONDITIONER;
2040: PetscInfo2(snes,"Reusing Jacobian/preconditioner because lag is %D and SNES iteration is %D\n",snes->lagjacobian,snes->iter);
2041: PetscObjectTypeCompare((PetscObject)*A,MATMFFD,&flag);
2042: if (flag) {
2043: MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
2044: MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
2045: }
2046: return(0);
2047: }
2049: *flg = DIFFERENT_NONZERO_PATTERN;
2050: PetscLogEventBegin(SNES_JacobianEval,snes,X,*A,*B);
2051: PetscStackPush("SNES user Jacobian function");
2052: (*sdm->computejacobian)(snes,X,A,B,flg,sdm->jacobianctx);
2053: PetscStackPop;
2054: PetscLogEventEnd(SNES_JacobianEval,snes,X,*A,*B);
2056: if (snes->lagpreconditioner == -2) {
2057: PetscInfo(snes,"Rebuilding preconditioner exactly once since lag is -2\n");
2058: snes->lagpreconditioner = -1;
2059: } else if (snes->lagpreconditioner == -1) {
2060: *flg = SAME_PRECONDITIONER;
2061: PetscInfo(snes,"Reusing preconditioner because lag is -1\n");
2062: } else if (snes->lagpreconditioner > 1 && snes->iter % snes->lagpreconditioner) {
2063: *flg = SAME_PRECONDITIONER;
2064: PetscInfo2(snes,"Reusing preconditioner because lag is %D and SNES iteration is %D\n",snes->lagpreconditioner,snes->iter);
2065: }
2067: /* make sure user returned a correct Jacobian and preconditioner */
2070: {
2071: PetscBool flag = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_operator = PETSC_FALSE;
2072: PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit",&flag,PETSC_NULL);
2073: PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw",&flag_draw,PETSC_NULL);
2074: PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw_contour",&flag_contour,PETSC_NULL);
2075: PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_operator",&flag_operator,PETSC_NULL);
2076: if (flag || flag_draw || flag_contour) {
2077: Mat Bexp_mine = PETSC_NULL,Bexp,FDexp;
2078: MatStructure mstruct;
2079: PetscViewer vdraw,vstdout;
2080: PetscBool flg;
2081: if (flag_operator) {
2082: MatComputeExplicitOperator(*A,&Bexp_mine);
2083: Bexp = Bexp_mine;
2084: } else {
2085: /* See if the preconditioning matrix can be viewed and added directly */
2086: PetscObjectTypeCompareAny((PetscObject)*B,&flg,MATSEQAIJ,MATMPIAIJ,MATSEQDENSE,MATMPIDENSE,MATSEQBAIJ,MATMPIBAIJ,MATSEQSBAIJ,MATMPIBAIJ,"");
2087: if (flg) Bexp = *B;
2088: else {
2089: /* If the "preconditioning" matrix is itself MATSHELL or some other type without direct support */
2090: MatComputeExplicitOperator(*B,&Bexp_mine);
2091: Bexp = Bexp_mine;
2092: }
2093: }
2094: MatConvert(Bexp,MATSAME,MAT_INITIAL_MATRIX,&FDexp);
2095: SNESDefaultComputeJacobian(snes,X,&FDexp,&FDexp,&mstruct,NULL);
2096: PetscViewerASCIIGetStdout(((PetscObject)snes)->comm,&vstdout);
2097: if (flag_draw || flag_contour) {
2098: PetscViewerDrawOpen(((PetscObject)snes)->comm,0,"Explicit Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);
2099: if (flag_contour) {PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);}
2100: } else vdraw = PETSC_NULL;
2101: PetscViewerASCIIPrintf(vstdout,"Explicit %s\n",flag_operator?"Jacobian":"preconditioning Jacobian");
2102: if (flag) {MatView(Bexp,vstdout);}
2103: if (vdraw) {MatView(Bexp,vdraw);}
2104: PetscViewerASCIIPrintf(vstdout,"Finite difference Jacobian\n");
2105: if (flag) {MatView(FDexp,vstdout);}
2106: if (vdraw) {MatView(FDexp,vdraw);}
2107: MatAYPX(FDexp,-1.0,Bexp,SAME_NONZERO_PATTERN);
2108: PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian\n");
2109: if (flag) {MatView(FDexp,vstdout);}
2110: if (vdraw) { /* Always use contour for the difference */
2111: PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);
2112: MatView(FDexp,vdraw);
2113: PetscViewerPopFormat(vdraw);
2114: }
2115: if (flag_contour) {PetscViewerPopFormat(vdraw);}
2116: PetscViewerDestroy(&vdraw);
2117: MatDestroy(&Bexp_mine);
2118: MatDestroy(&FDexp);
2119: }
2120: }
2121: {
2122: PetscBool flag = PETSC_FALSE,flag_display = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_threshold = PETSC_FALSE;
2123: PetscReal threshold_atol = PETSC_SQRT_MACHINE_EPSILON,threshold_rtol = 10*PETSC_SQRT_MACHINE_EPSILON;
2124: PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring",&flag,PETSC_NULL);
2125: PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_display",&flag_display,PETSC_NULL);
2126: PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw",&flag_draw,PETSC_NULL);
2127: PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw_contour",&flag_contour,PETSC_NULL);
2128: PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold",&flag_threshold,PETSC_NULL);
2129: PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_rtol",&threshold_rtol,PETSC_NULL);
2130: PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_atol",&threshold_atol,PETSC_NULL);
2131: if (flag || flag_display || flag_draw || flag_contour || flag_threshold) {
2132: Mat Bfd;
2133: MatStructure mstruct;
2134: PetscViewer vdraw,vstdout;
2135: ISColoring iscoloring;
2136: MatFDColoring matfdcoloring;
2137: PetscErrorCode (*func)(SNES,Vec,Vec,void*);
2138: void *funcctx;
2139: PetscReal norm1,norm2,normmax;
2141: MatDuplicate(*B,MAT_DO_NOT_COPY_VALUES,&Bfd);
2142: MatGetColoring(Bfd,MATCOLORINGSL,&iscoloring);
2143: MatFDColoringCreate(Bfd,iscoloring,&matfdcoloring);
2144: ISColoringDestroy(&iscoloring);
2146: /* This method of getting the function is currently unreliable since it doesn't work for DM local functions. */
2147: SNESGetFunction(snes,PETSC_NULL,&func,&funcctx);
2148: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode(*)(void))func,funcctx);
2149: PetscObjectSetOptionsPrefix((PetscObject)matfdcoloring,((PetscObject)snes)->prefix);
2150: PetscObjectAppendOptionsPrefix((PetscObject)matfdcoloring,"coloring_");
2151: MatFDColoringSetFromOptions(matfdcoloring);
2152: MatFDColoringApply(Bfd,matfdcoloring,X,&mstruct,snes);
2153: MatFDColoringDestroy(&matfdcoloring);
2155: PetscViewerASCIIGetStdout(((PetscObject)snes)->comm,&vstdout);
2156: if (flag_draw || flag_contour) {
2157: PetscViewerDrawOpen(((PetscObject)snes)->comm,0,"Colored Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);
2158: if (flag_contour) {PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);}
2159: } else vdraw = PETSC_NULL;
2160: PetscViewerASCIIPrintf(vstdout,"Explicit preconditioning Jacobian\n");
2161: if (flag_display) {MatView(*B,vstdout);}
2162: if (vdraw) {MatView(*B,vdraw);}
2163: PetscViewerASCIIPrintf(vstdout,"Colored Finite difference Jacobian\n");
2164: if (flag_display) {MatView(Bfd,vstdout);}
2165: if (vdraw) {MatView(Bfd,vdraw);}
2166: MatAYPX(Bfd,-1.0,*B,SAME_NONZERO_PATTERN);
2167: MatNorm(Bfd,NORM_1,&norm1);
2168: MatNorm(Bfd,NORM_FROBENIUS,&norm2);
2169: MatNorm(Bfd,NORM_MAX,&normmax);
2170: PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian, norm1=%G normFrob=%G normmax=%G\n",norm1,norm2,normmax);
2171: if (flag_display) {MatView(Bfd,vstdout);}
2172: if (vdraw) { /* Always use contour for the difference */
2173: PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);
2174: MatView(Bfd,vdraw);
2175: PetscViewerPopFormat(vdraw);
2176: }
2177: if (flag_contour) {PetscViewerPopFormat(vdraw);}
2179: if (flag_threshold) {
2180: PetscInt bs,rstart,rend,i;
2181: MatGetBlockSize(*B,&bs);
2182: MatGetOwnershipRange(*B,&rstart,&rend);
2183: for (i=rstart; i<rend; i++) {
2184: const PetscScalar *ba,*ca;
2185: const PetscInt *bj,*cj;
2186: PetscInt bn,cn,j,maxentrycol = -1,maxdiffcol = -1,maxrdiffcol = -1;
2187: PetscReal maxentry = 0,maxdiff = 0,maxrdiff = 0;
2188: MatGetRow(*B,i,&bn,&bj,&ba);
2189: MatGetRow(Bfd,i,&cn,&cj,&ca);
2190: if (bn != cn) SETERRQ(((PetscObject)*A)->comm,PETSC_ERR_PLIB,"Unexpected different nonzero pattern in -snes_compare_coloring_threshold");
2191: for (j=0; j<bn; j++) {
2192: PetscReal rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j]));
2193: if (PetscAbsScalar(ba[j]) > PetscAbs(maxentry)) {
2194: maxentrycol = bj[j];
2195: maxentry = PetscRealPart(ba[j]);
2196: }
2197: if (PetscAbsScalar(ca[j]) > PetscAbs(maxdiff)) {
2198: maxdiffcol = bj[j];
2199: maxdiff = PetscRealPart(ca[j]);
2200: }
2201: if (rdiff > maxrdiff) {
2202: maxrdiffcol = bj[j];
2203: maxrdiff = rdiff;
2204: }
2205: }
2206: if (maxrdiff > 1) {
2207: PetscViewerASCIIPrintf(vstdout,"row %D (maxentry=%G at %D, maxdiff=%G at %D, maxrdiff=%G at %D):",i,maxentry,maxentrycol,maxdiff,maxdiffcol,maxrdiff,maxrdiffcol);
2208: for (j=0; j<bn; j++) {
2209: PetscReal rdiff;
2210: rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j]));
2211: if (rdiff > 1) {
2212: PetscViewerASCIIPrintf(vstdout," (%D,%G:%G)",bj[j],PetscRealPart(ba[j]),PetscRealPart(ca[j]));
2213: }
2214: }
2215: PetscViewerASCIIPrintf(vstdout,"\n",i,maxentry,maxdiff,maxrdiff);
2216: }
2217: MatRestoreRow(*B,i,&bn,&bj,&ba);
2218: MatRestoreRow(Bfd,i,&cn,&cj,&ca);
2219: }
2220: }
2221: PetscViewerDestroy(&vdraw);
2222: MatDestroy(&Bfd);
2223: }
2224: }
2225: return(0);
2226: }
2230: /*@C
2231: SNESSetJacobian - Sets the function to compute Jacobian as well as the
2232: location to store the matrix.
2234: Logically Collective on SNES and Mat
2236: Input Parameters:
2237: + snes - the SNES context
2238: . A - Jacobian matrix
2239: . B - preconditioner matrix (usually same as the Jacobian)
2240: . func - Jacobian evaluation routine (if PETSC_NULL then SNES retains any previously set value)
2241: - ctx - [optional] user-defined context for private data for the
2242: Jacobian evaluation routine (may be PETSC_NULL) (if PETSC_NULL then SNES retains any previously set value)
2244: Calling sequence of func:
2245: $ func (SNES snes,Vec x,Mat *A,Mat *B,int *flag,void *ctx);
2247: + x - input vector
2248: . A - Jacobian matrix
2249: . B - preconditioner matrix, usually the same as A
2250: . flag - flag indicating information about the preconditioner matrix
2251: structure (same as flag in KSPSetOperators()), one of SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER
2252: - ctx - [optional] user-defined Jacobian context
2254: Notes:
2255: See KSPSetOperators() for important information about setting the flag
2256: output parameter in the routine func(). Be sure to read this information!
2258: The routine func() takes Mat * as the matrix arguments rather than Mat.
2259: This allows the Jacobian evaluation routine to replace A and/or B with a
2260: completely new new matrix structure (not just different matrix elements)
2261: when appropriate, for instance, if the nonzero structure is changing
2262: throughout the global iterations.
2264: If the A matrix and B matrix are different you must call MatAssemblyBegin/End() on
2265: each matrix.
2267: If using SNESDefaultComputeJacobianColor() to assemble a Jacobian, the ctx argument
2268: must be a MatFDColoring.
2270: Other defect-correction schemes can be used by computing a different matrix in place of the Jacobian. One common
2271: example is to use the "Picard linearization" which only differentiates through the highest order parts of each term.
2273: Level: beginner
2275: .keywords: SNES, nonlinear, set, Jacobian, matrix
2277: .seealso: KSPSetOperators(), SNESSetFunction(), MatMFFDComputeJacobian(), SNESDefaultComputeJacobianColor(), MatStructure
2278: @*/
2279: PetscErrorCode SNESSetJacobian(SNES snes,Mat A,Mat B,PetscErrorCode (*func)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void *ctx)
2280: {
2282: DM dm;
2290: SNESGetDM(snes,&dm);
2291: DMSNESSetJacobian(dm,func,ctx);
2292: if (A) {
2293: PetscObjectReference((PetscObject)A);
2294: MatDestroy(&snes->jacobian);
2295: snes->jacobian = A;
2296: }
2297: if (B) {
2298: PetscObjectReference((PetscObject)B);
2299: MatDestroy(&snes->jacobian_pre);
2300: snes->jacobian_pre = B;
2301: }
2302: return(0);
2303: }
2307: /*@C
2308: SNESGetJacobian - Returns the Jacobian matrix and optionally the user
2309: provided context for evaluating the Jacobian.
2311: Not Collective, but Mat object will be parallel if SNES object is
2313: Input Parameter:
2314: . snes - the nonlinear solver context
2316: Output Parameters:
2317: + A - location to stash Jacobian matrix (or PETSC_NULL)
2318: . B - location to stash preconditioner matrix (or PETSC_NULL)
2319: . func - location to put Jacobian function (or PETSC_NULL)
2320: - ctx - location to stash Jacobian ctx (or PETSC_NULL)
2322: Level: advanced
2324: .seealso: SNESSetJacobian(), SNESComputeJacobian()
2325: @*/
2326: PetscErrorCode SNESGetJacobian(SNES snes,Mat *A,Mat *B,PetscErrorCode (**func)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void **ctx)
2327: {
2329: DM dm;
2330: SNESDM sdm;
2334: if (A) *A = snes->jacobian;
2335: if (B) *B = snes->jacobian_pre;
2336: SNESGetDM(snes,&dm);
2337: DMSNESGetContext(dm,&sdm);
2338: if (func) *func = sdm->computejacobian;
2339: if (ctx) *ctx = sdm->jacobianctx;
2340: return(0);
2341: }
2343: /* ----- Routines to initialize and destroy a nonlinear solver ---- */
2347: /*@
2348: SNESSetUp - Sets up the internal data structures for the later use
2349: of a nonlinear solver.
2351: Collective on SNES
2353: Input Parameters:
2354: . snes - the SNES context
2356: Notes:
2357: For basic use of the SNES solvers the user need not explicitly call
2358: SNESSetUp(), since these actions will automatically occur during
2359: the call to SNESSolve(). However, if one wishes to control this
2360: phase separately, SNESSetUp() should be called after SNESCreate()
2361: and optional routines of the form SNESSetXXX(), but before SNESSolve().
2363: Level: advanced
2365: .keywords: SNES, nonlinear, setup
2367: .seealso: SNESCreate(), SNESSolve(), SNESDestroy()
2368: @*/
2369: PetscErrorCode SNESSetUp(SNES snes)
2370: {
2372: DM dm;
2373: SNESDM sdm;
2374: SNESLineSearch linesearch;
2375: SNESLineSearch pclinesearch;
2376: void *lsprectx,*lspostctx;
2377: SNESLineSearchPreCheckFunc lsprefunc;
2378: SNESLineSearchPostCheckFunc lspostfunc;
2379: PetscErrorCode (*func)(SNES,Vec,Vec,void*);
2380: Vec f,fpc;
2381: void *funcctx;
2382: PetscErrorCode (*jac)(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
2383: void *jacctx;
2384: Mat A,B;
2388: if (snes->setupcalled) return(0);
2390: if (!((PetscObject)snes)->type_name) {
2391: SNESSetType(snes,SNESLS);
2392: }
2394: SNESGetFunction(snes,&snes->vec_func,PETSC_NULL,PETSC_NULL);
2395: if (snes->vec_func == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be function vector");
2396: if (snes->vec_rhs == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be right hand side vector");
2398: if (!snes->vec_sol_update /* && snes->vec_sol */) {
2399: VecDuplicate(snes->vec_sol,&snes->vec_sol_update);
2400: PetscLogObjectParent(snes,snes->vec_sol_update);
2401: }
2403: SNESGetDM(snes,&dm);
2404: DMSNESSetUpLegacy(dm); /* To be removed when function routines are taken out of the DM package */
2405: DMSNESGetContext(dm,&sdm);
2406: if (!sdm->computefunction) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must provide a residual function with SNESSetFunction(), DMSNESSetFunction(), DMDASNESSetFunctionLocal(), etc");
2407: if (!snes->vec_func) {
2408: DMCreateGlobalVector(dm,&snes->vec_func);
2409: }
2411: if (!snes->ksp) {SNESGetKSP(snes, &snes->ksp);}
2413: if (!snes->linesearch) {SNESGetSNESLineSearch(snes, &snes->linesearch);}
2415: if (snes->ops->usercompute && !snes->user) {
2416: (*snes->ops->usercompute)(snes,(void**)&snes->user);
2417: }
2419: if (snes->pc) {
2420: /* copy the DM over */
2421: SNESGetDM(snes,&dm);
2422: SNESSetDM(snes->pc,dm);
2424: /* copy the legacy SNES context not related to the DM over*/
2425: SNESGetFunction(snes,&f,&func,&funcctx);
2426: VecDuplicate(f,&fpc);
2427: SNESSetFunction(snes->pc,fpc,func,funcctx);
2428: SNESGetJacobian(snes,&A,&B,&jac,&jacctx);
2429: SNESSetJacobian(snes->pc,A,B,jac,jacctx);
2430: VecDestroy(&fpc);
2432: /* copy the function pointers over */
2433: PetscObjectCopyFortranFunctionPointers((PetscObject)snes,(PetscObject)snes->pc);
2435: /* default to 1 iteration */
2436: SNESSetTolerances(snes->pc, 0.0, 0.0, 0.0, 1, snes->pc->max_funcs);
2437: SNESSetNormType(snes->pc, SNES_NORM_FINAL_ONLY);
2438: SNESSetFromOptions(snes->pc);
2440: /* copy the line search context over */
2441: SNESGetSNESLineSearch(snes,&linesearch);
2442: SNESGetSNESLineSearch(snes->pc,&pclinesearch);
2443: SNESLineSearchGetPreCheck(linesearch,&lsprefunc,&lsprectx);
2444: SNESLineSearchGetPostCheck(linesearch,&lspostfunc,&lspostctx);
2445: SNESLineSearchSetPreCheck(pclinesearch,lsprefunc,lsprectx);
2446: SNESLineSearchSetPostCheck(pclinesearch,lspostfunc,lspostctx);
2447: PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)pclinesearch);
2448: }
2450: if (snes->ops->setup) {
2451: (*snes->ops->setup)(snes);
2452: }
2454: snes->setupcalled = PETSC_TRUE;
2455: return(0);
2456: }
2460: /*@
2461: SNESReset - Resets a SNES context to the snessetupcalled = 0 state and removes any allocated Vecs and Mats
2463: Collective on SNES
2465: Input Parameter:
2466: . snes - iterative context obtained from SNESCreate()
2468: Level: intermediate
2470: Notes: Also calls the application context destroy routine set with SNESSetComputeApplicationContext()
2472: .keywords: SNES, destroy
2474: .seealso: SNESCreate(), SNESSetUp(), SNESSolve()
2475: @*/
2476: PetscErrorCode SNESReset(SNES snes)
2477: {
2482: if (snes->ops->userdestroy && snes->user) {
2483: (*snes->ops->userdestroy)((void**)&snes->user);
2484: snes->user = PETSC_NULL;
2485: }
2486: if (snes->pc) {
2487: SNESReset(snes->pc);
2488: }
2490: if (snes->ops->reset) {
2491: (*snes->ops->reset)(snes);
2492: }
2493: if (snes->ksp) {
2494: KSPReset(snes->ksp);
2495: }
2497: if (snes->linesearch) {
2498: SNESLineSearchReset(snes->linesearch);
2499: }
2501: VecDestroy(&snes->vec_rhs);
2502: VecDestroy(&snes->vec_sol);
2503: VecDestroy(&snes->vec_sol_update);
2504: VecDestroy(&snes->vec_func);
2505: MatDestroy(&snes->jacobian);
2506: MatDestroy(&snes->jacobian_pre);
2507: VecDestroyVecs(snes->nwork,&snes->work);
2508: VecDestroyVecs(snes->nvwork,&snes->vwork);
2509: snes->nwork = snes->nvwork = 0;
2510: snes->setupcalled = PETSC_FALSE;
2511: return(0);
2512: }
2516: /*@
2517: SNESDestroy - Destroys the nonlinear solver context that was created
2518: with SNESCreate().
2520: Collective on SNES
2522: Input Parameter:
2523: . snes - the SNES context
2525: Level: beginner
2527: .keywords: SNES, nonlinear, destroy
2529: .seealso: SNESCreate(), SNESSolve()
2530: @*/
2531: PetscErrorCode SNESDestroy(SNES *snes)
2532: {
2536: if (!*snes) return(0);
2538: if (--((PetscObject)(*snes))->refct > 0) {*snes = 0; return(0);}
2540: SNESReset((*snes));
2541: SNESDestroy(&(*snes)->pc);
2543: /* if memory was published with AMS then destroy it */
2544: PetscObjectDepublish((*snes));
2545: if ((*snes)->ops->destroy) {(*((*snes))->ops->destroy)((*snes));}
2547: DMDestroy(&(*snes)->dm);
2548: KSPDestroy(&(*snes)->ksp);
2549: SNESLineSearchDestroy(&(*snes)->linesearch);
2551: PetscFree((*snes)->kspconvctx);
2552: if ((*snes)->ops->convergeddestroy) {
2553: (*(*snes)->ops->convergeddestroy)((*snes)->cnvP);
2554: }
2555: if ((*snes)->conv_malloc) {
2556: PetscFree((*snes)->conv_hist);
2557: PetscFree((*snes)->conv_hist_its);
2558: }
2559: SNESMonitorCancel((*snes));
2560: PetscHeaderDestroy(snes);
2561: return(0);
2562: }
2564: /* ----------- Routines to set solver parameters ---------- */
2568: /*@
2569: SNESSetLagPreconditioner - Determines when the preconditioner is rebuilt in the nonlinear solve.
2571: Logically Collective on SNES
2573: Input Parameters:
2574: + snes - the SNES context
2575: - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time
2576: the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that
2578: Options Database Keys:
2579: . -snes_lag_preconditioner <lag>
2581: Notes:
2582: The default is 1
2583: The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1
2584: If -1 is used before the very first nonlinear solve the preconditioner is still built because there is no previous preconditioner to use
2586: Level: intermediate
2588: .keywords: SNES, nonlinear, set, convergence, tolerances
2590: .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian()
2592: @*/
2593: PetscErrorCode SNESSetLagPreconditioner(SNES snes,PetscInt lag)
2594: {
2597: if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater");
2598: if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0");
2600: snes->lagpreconditioner = lag;
2601: return(0);
2602: }
2606: /*@
2607: SNESSetGridSequence - sets the number of steps of grid sequencing that SNES does
2609: Logically Collective on SNES
2611: Input Parameters:
2612: + snes - the SNES context
2613: - steps - the number of refinements to do, defaults to 0
2615: Options Database Keys:
2616: . -snes_grid_sequence <steps>
2618: Level: intermediate
2620: Notes:
2621: Use SNESGetSolution() to extract the fine grid solution after grid sequencing.
2623: .keywords: SNES, nonlinear, set, convergence, tolerances
2625: .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian()
2627: @*/
2628: PetscErrorCode SNESSetGridSequence(SNES snes,PetscInt steps)
2629: {
2633: snes->gridsequence = steps;
2634: return(0);
2635: }
2639: /*@
2640: SNESGetLagPreconditioner - Indicates how often the preconditioner is rebuilt
2642: Not Collective
2644: Input Parameter:
2645: . snes - the SNES context
2646:
2647: Output Parameter:
2648: . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time
2649: the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that
2651: Options Database Keys:
2652: . -snes_lag_preconditioner <lag>
2654: Notes:
2655: The default is 1
2656: The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1
2658: Level: intermediate
2660: .keywords: SNES, nonlinear, set, convergence, tolerances
2662: .seealso: SNESSetTrustRegionTolerance(), SNESSetLagPreconditioner()
2664: @*/
2665: PetscErrorCode SNESGetLagPreconditioner(SNES snes,PetscInt *lag)
2666: {
2669: *lag = snes->lagpreconditioner;
2670: return(0);
2671: }
2675: /*@
2676: SNESSetLagJacobian - Determines when the Jacobian is rebuilt in the nonlinear solve. See SNESSetLagPreconditioner() for determining how
2677: often the preconditioner is rebuilt.
2679: Logically Collective on SNES
2681: Input Parameters:
2682: + snes - the SNES context
2683: - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time
2684: the Jacobian is built etc. -2 means rebuild at next chance but then never again
2686: Options Database Keys:
2687: . -snes_lag_jacobian <lag>
2689: Notes:
2690: The default is 1
2691: The Jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1
2692: 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
2693: at the next Newton step but never again (unless it is reset to another value)
2695: Level: intermediate
2697: .keywords: SNES, nonlinear, set, convergence, tolerances
2699: .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagPreconditioner(), SNESGetLagJacobian()
2701: @*/
2702: PetscErrorCode SNESSetLagJacobian(SNES snes,PetscInt lag)
2703: {
2706: if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater");
2707: if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0");
2709: snes->lagjacobian = lag;
2710: return(0);
2711: }
2715: /*@
2716: SNESGetLagJacobian - Indicates how often the Jacobian is rebuilt. See SNESGetLagPreconditioner() to determine when the preconditioner is rebuilt
2718: Not Collective
2720: Input Parameter:
2721: . snes - the SNES context
2722:
2723: Output Parameter:
2724: . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time
2725: the Jacobian is built etc.
2727: Options Database Keys:
2728: . -snes_lag_jacobian <lag>
2730: Notes:
2731: The default is 1
2732: The jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1
2734: Level: intermediate
2736: .keywords: SNES, nonlinear, set, convergence, tolerances
2738: .seealso: SNESSetTrustRegionTolerance(), SNESSetLagJacobian(), SNESSetLagPreconditioner(), SNESGetLagPreconditioner()
2740: @*/
2741: PetscErrorCode SNESGetLagJacobian(SNES snes,PetscInt *lag)
2742: {
2745: *lag = snes->lagjacobian;
2746: return(0);
2747: }
2751: /*@
2752: SNESSetTolerances - Sets various parameters used in convergence tests.
2754: Logically Collective on SNES
2756: Input Parameters:
2757: + snes - the SNES context
2758: . abstol - absolute convergence tolerance
2759: . rtol - relative convergence tolerance
2760: . stol - convergence tolerance in terms of the norm
2761: of the change in the solution between steps
2762: . maxit - maximum number of iterations
2763: - maxf - maximum number of function evaluations
2765: Options Database Keys:
2766: + -snes_atol <abstol> - Sets abstol
2767: . -snes_rtol <rtol> - Sets rtol
2768: . -snes_stol <stol> - Sets stol
2769: . -snes_max_it <maxit> - Sets maxit
2770: - -snes_max_funcs <maxf> - Sets maxf
2772: Notes:
2773: The default maximum number of iterations is 50.
2774: The default maximum number of function evaluations is 1000.
2776: Level: intermediate
2778: .keywords: SNES, nonlinear, set, convergence, tolerances
2780: .seealso: SNESSetTrustRegionTolerance()
2781: @*/
2782: PetscErrorCode SNESSetTolerances(SNES snes,PetscReal abstol,PetscReal rtol,PetscReal stol,PetscInt maxit,PetscInt maxf)
2783: {
2792: if (abstol != PETSC_DEFAULT) {
2793: if (abstol < 0.0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %G must be non-negative",abstol);
2794: snes->abstol = abstol;
2795: }
2796: if (rtol != PETSC_DEFAULT) {
2797: 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);
2798: snes->rtol = rtol;
2799: }
2800: if (stol != PETSC_DEFAULT) {
2801: if (stol < 0.0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Step tolerance %G must be non-negative",stol);
2802: snes->stol = stol;
2803: }
2804: if (maxit != PETSC_DEFAULT) {
2805: if (maxit < 0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxit);
2806: snes->max_its = maxit;
2807: }
2808: if (maxf != PETSC_DEFAULT) {
2809: if (maxf < 0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of function evaluations %D must be non-negative",maxf);
2810: snes->max_funcs = maxf;
2811: }
2812: snes->tolerancesset = PETSC_TRUE;
2813: return(0);
2814: }
2818: /*@
2819: SNESGetTolerances - Gets various parameters used in convergence tests.
2821: Not Collective
2823: Input Parameters:
2824: + snes - the SNES context
2825: . atol - absolute convergence tolerance
2826: . rtol - relative convergence tolerance
2827: . stol - convergence tolerance in terms of the norm
2828: of the change in the solution between steps
2829: . maxit - maximum number of iterations
2830: - maxf - maximum number of function evaluations
2832: Notes:
2833: The user can specify PETSC_NULL for any parameter that is not needed.
2835: Level: intermediate
2837: .keywords: SNES, nonlinear, get, convergence, tolerances
2839: .seealso: SNESSetTolerances()
2840: @*/
2841: PetscErrorCode SNESGetTolerances(SNES snes,PetscReal *atol,PetscReal *rtol,PetscReal *stol,PetscInt *maxit,PetscInt *maxf)
2842: {
2845: if (atol) *atol = snes->abstol;
2846: if (rtol) *rtol = snes->rtol;
2847: if (stol) *stol = snes->stol;
2848: if (maxit) *maxit = snes->max_its;
2849: if (maxf) *maxf = snes->max_funcs;
2850: return(0);
2851: }
2855: /*@
2856: SNESSetTrustRegionTolerance - Sets the trust region parameter tolerance.
2858: Logically Collective on SNES
2860: Input Parameters:
2861: + snes - the SNES context
2862: - tol - tolerance
2863:
2864: Options Database Key:
2865: . -snes_trtol <tol> - Sets tol
2867: Level: intermediate
2869: .keywords: SNES, nonlinear, set, trust region, tolerance
2871: .seealso: SNESSetTolerances()
2872: @*/
2873: PetscErrorCode SNESSetTrustRegionTolerance(SNES snes,PetscReal tol)
2874: {
2878: snes->deltatol = tol;
2879: return(0);
2880: }
2882: /*
2883: Duplicate the lg monitors for SNES from KSP; for some reason with
2884: dynamic libraries things don't work under Sun4 if we just use
2885: macros instead of functions
2886: */
2889: PetscErrorCode SNESMonitorLG(SNES snes,PetscInt it,PetscReal norm,void *ctx)
2890: {
2895: KSPMonitorLG((KSP)snes,it,norm,ctx);
2896: return(0);
2897: }
2901: PetscErrorCode SNESMonitorLGCreate(const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *draw)
2902: {
2906: KSPMonitorLGCreate(host,label,x,y,m,n,draw);
2907: return(0);
2908: }
2912: PetscErrorCode SNESMonitorLGDestroy(PetscDrawLG *draw)
2913: {
2917: KSPMonitorLGDestroy(draw);
2918: return(0);
2919: }
2921: extern PetscErrorCode SNESMonitorRange_Private(SNES,PetscInt,PetscReal*);
2924: PetscErrorCode SNESMonitorLGRange(SNES snes,PetscInt n,PetscReal rnorm,void *monctx)
2925: {
2926: PetscDrawLG lg;
2927: PetscErrorCode ierr;
2928: PetscReal x,y,per;
2929: PetscViewer v = (PetscViewer)monctx;
2930: static PetscReal prev; /* should be in the context */
2931: PetscDraw draw;
2933: if (!monctx) {
2934: MPI_Comm comm;
2936: PetscObjectGetComm((PetscObject)snes,&comm);
2937: v = PETSC_VIEWER_DRAW_(comm);
2938: }
2939: PetscViewerDrawGetDrawLG(v,0,&lg);
2940: if (!n) {PetscDrawLGReset(lg);}
2941: PetscDrawLGGetDraw(lg,&draw);
2942: PetscDrawSetTitle(draw,"Residual norm");
2943: x = (PetscReal) n;
2944: if (rnorm > 0.0) y = log10(rnorm); else y = -15.0;
2945: PetscDrawLGAddPoint(lg,&x,&y);
2946: if (n < 20 || !(n % 5)) {
2947: PetscDrawLGDraw(lg);
2948: }
2950: PetscViewerDrawGetDrawLG(v,1,&lg);
2951: if (!n) {PetscDrawLGReset(lg);}
2952: PetscDrawLGGetDraw(lg,&draw);
2953: PetscDrawSetTitle(draw,"% elemts > .2*max elemt");
2954: SNESMonitorRange_Private(snes,n,&per);
2955: x = (PetscReal) n;
2956: y = 100.0*per;
2957: PetscDrawLGAddPoint(lg,&x,&y);
2958: if (n < 20 || !(n % 5)) {
2959: PetscDrawLGDraw(lg);
2960: }
2962: PetscViewerDrawGetDrawLG(v,2,&lg);
2963: if (!n) {prev = rnorm;PetscDrawLGReset(lg);}
2964: PetscDrawLGGetDraw(lg,&draw);
2965: PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm");
2966: x = (PetscReal) n;
2967: y = (prev - rnorm)/prev;
2968: PetscDrawLGAddPoint(lg,&x,&y);
2969: if (n < 20 || !(n % 5)) {
2970: PetscDrawLGDraw(lg);
2971: }
2973: PetscViewerDrawGetDrawLG(v,3,&lg);
2974: if (!n) {PetscDrawLGReset(lg);}
2975: PetscDrawLGGetDraw(lg,&draw);
2976: PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm*(% > .2 max)");
2977: x = (PetscReal) n;
2978: y = (prev - rnorm)/(prev*per);
2979: if (n > 2) { /*skip initial crazy value */
2980: PetscDrawLGAddPoint(lg,&x,&y);
2981: }
2982: if (n < 20 || !(n % 5)) {
2983: PetscDrawLGDraw(lg);
2984: }
2985: prev = rnorm;
2986: return(0);
2987: }
2991: PetscErrorCode SNESMonitorLGRangeCreate(const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *draw)
2992: {
2996: KSPMonitorLGCreate(host,label,x,y,m,n,draw);
2997: return(0);
2998: }
3002: PetscErrorCode SNESMonitorLGRangeDestroy(PetscDrawLG *draw)
3003: {
3007: KSPMonitorLGDestroy(draw);
3008: return(0);
3009: }
3013: /*@
3014: SNESMonitor - runs the user provided monitor routines, if they exist
3016: Collective on SNES
3018: Input Parameters:
3019: + snes - nonlinear solver context obtained from SNESCreate()
3020: . iter - iteration number
3021: - rnorm - relative norm of the residual
3023: Notes:
3024: This routine is called by the SNES implementations.
3025: It does not typically need to be called by the user.
3027: Level: developer
3029: .seealso: SNESMonitorSet()
3030: @*/
3031: PetscErrorCode SNESMonitor(SNES snes,PetscInt iter,PetscReal rnorm)
3032: {
3034: PetscInt i,n = snes->numbermonitors;
3037: for (i=0; i<n; i++) {
3038: (*snes->monitor[i])(snes,iter,rnorm,snes->monitorcontext[i]);
3039: }
3040: return(0);
3041: }
3043: /* ------------ Routines to set performance monitoring options ----------- */
3047: /*@C
3048: SNESMonitorSet - Sets an ADDITIONAL function that is to be used at every
3049: iteration of the nonlinear solver to display the iteration's
3050: progress.
3052: Logically Collective on SNES
3054: Input Parameters:
3055: + snes - the SNES context
3056: . func - monitoring routine
3057: . mctx - [optional] user-defined context for private data for the
3058: monitor routine (use PETSC_NULL if no context is desired)
3059: - monitordestroy - [optional] routine that frees monitor context
3060: (may be PETSC_NULL)
3062: Calling sequence of func:
3063: $ int func(SNES snes,PetscInt its, PetscReal norm,void *mctx)
3065: + snes - the SNES context
3066: . its - iteration number
3067: . norm - 2-norm function value (may be estimated)
3068: - mctx - [optional] monitoring context
3070: Options Database Keys:
3071: + -snes_monitor - sets SNESMonitorDefault()
3072: . -snes_monitor_draw - sets line graph monitor,
3073: uses SNESMonitorLGCreate()
3074: - -snes_monitor_cancel - cancels all monitors that have
3075: been hardwired into a code by
3076: calls to SNESMonitorSet(), but
3077: does not cancel those set via
3078: the options database.
3080: Notes:
3081: Several different monitoring routines may be set by calling
3082: SNESMonitorSet() multiple times; all will be called in the
3083: order in which they were set.
3085: Fortran notes: Only a single monitor function can be set for each SNES object
3087: Level: intermediate
3089: .keywords: SNES, nonlinear, set, monitor
3091: .seealso: SNESMonitorDefault(), SNESMonitorCancel()
3092: @*/
3093: PetscErrorCode SNESMonitorSet(SNES snes,PetscErrorCode (*monitor)(SNES,PetscInt,PetscReal,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**))
3094: {
3095: PetscInt i;
3100: if (snes->numbermonitors >= MAXSNESMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set");
3101: for (i=0; i<snes->numbermonitors;i++) {
3102: if (monitor == snes->monitor[i] && monitordestroy == snes->monitordestroy[i] && mctx == snes->monitorcontext[i]) {
3103: if (monitordestroy) {
3104: (*monitordestroy)(&mctx);
3105: }
3106: return(0);
3107: }
3108: }
3109: snes->monitor[snes->numbermonitors] = monitor;
3110: snes->monitordestroy[snes->numbermonitors] = monitordestroy;
3111: snes->monitorcontext[snes->numbermonitors++] = (void*)mctx;
3112: return(0);
3113: }
3117: /*@C
3118: SNESMonitorCancel - Clears all the monitor functions for a SNES object.
3120: Logically Collective on SNES
3122: Input Parameters:
3123: . snes - the SNES context
3125: Options Database Key:
3126: . -snes_monitor_cancel - cancels all monitors that have been hardwired
3127: into a code by calls to SNESMonitorSet(), but does not cancel those
3128: set via the options database
3130: Notes:
3131: There is no way to clear one specific monitor from a SNES object.
3133: Level: intermediate
3135: .keywords: SNES, nonlinear, set, monitor
3137: .seealso: SNESMonitorDefault(), SNESMonitorSet()
3138: @*/
3139: PetscErrorCode SNESMonitorCancel(SNES snes)
3140: {
3142: PetscInt i;
3146: for (i=0; i<snes->numbermonitors; i++) {
3147: if (snes->monitordestroy[i]) {
3148: (*snes->monitordestroy[i])(&snes->monitorcontext[i]);
3149: }
3150: }
3151: snes->numbermonitors = 0;
3152: return(0);
3153: }
3157: /*@C
3158: SNESSetConvergenceTest - Sets the function that is to be used
3159: to test for convergence of the nonlinear iterative solution.
3161: Logically Collective on SNES
3163: Input Parameters:
3164: + snes - the SNES context
3165: . func - routine to test for convergence
3166: . cctx - [optional] context for private data for the convergence routine (may be PETSC_NULL)
3167: - destroy - [optional] destructor for the context (may be PETSC_NULL; PETSC_NULL_FUNCTION in Fortran)
3169: Calling sequence of func:
3170: $ PetscErrorCode func (SNES snes,PetscInt it,PetscReal xnorm,PetscReal gnorm,PetscReal f,SNESConvergedReason *reason,void *cctx)
3172: + snes - the SNES context
3173: . it - current iteration (0 is the first and is before any Newton step)
3174: . cctx - [optional] convergence context
3175: . reason - reason for convergence/divergence
3176: . xnorm - 2-norm of current iterate
3177: . gnorm - 2-norm of current step
3178: - f - 2-norm of function
3180: Level: advanced
3182: .keywords: SNES, nonlinear, set, convergence, test
3184: .seealso: SNESDefaultConverged(), SNESSkipConverged()
3185: @*/
3186: PetscErrorCode SNESSetConvergenceTest(SNES snes,PetscErrorCode (*func)(SNES,PetscInt,PetscReal,PetscReal,PetscReal,SNESConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*))
3187: {
3192: if (!func) func = SNESSkipConverged;
3193: if (snes->ops->convergeddestroy) {
3194: (*snes->ops->convergeddestroy)(snes->cnvP);
3195: }
3196: snes->ops->converged = func;
3197: snes->ops->convergeddestroy = destroy;
3198: snes->cnvP = cctx;
3199: return(0);
3200: }
3204: /*@
3205: SNESGetConvergedReason - Gets the reason the SNES iteration was stopped.
3207: Not Collective
3209: Input Parameter:
3210: . snes - the SNES context
3212: Output Parameter:
3213: . reason - negative value indicates diverged, positive value converged, see SNESConvergedReason or the
3214: manual pages for the individual convergence tests for complete lists
3216: Level: intermediate
3218: Notes: Can only be called after the call the SNESSolve() is complete.
3220: .keywords: SNES, nonlinear, set, convergence, test
3222: .seealso: SNESSetConvergenceTest(), SNESConvergedReason
3223: @*/
3224: PetscErrorCode SNESGetConvergedReason(SNES snes,SNESConvergedReason *reason)
3225: {
3229: *reason = snes->reason;
3230: return(0);
3231: }
3235: /*@
3236: SNESSetConvergenceHistory - Sets the array used to hold the convergence history.
3238: Logically Collective on SNES
3240: Input Parameters:
3241: + snes - iterative context obtained from SNESCreate()
3242: . a - array to hold history, this array will contain the function norms computed at each step
3243: . its - integer array holds the number of linear iterations for each solve.
3244: . na - size of a and its
3245: - reset - PETSC_TRUE indicates each new nonlinear solve resets the history counter to zero,
3246: else it continues storing new values for new nonlinear solves after the old ones
3248: Notes:
3249: If 'a' and 'its' are PETSC_NULL then space is allocated for the history. If 'na' PETSC_DECIDE or PETSC_DEFAULT then a
3250: default array of length 10000 is allocated.
3252: This routine is useful, e.g., when running a code for purposes
3253: of accurate performance monitoring, when no I/O should be done
3254: during the section of code that is being timed.
3256: Level: intermediate
3258: .keywords: SNES, set, convergence, history
3260: .seealso: SNESGetConvergenceHistory()
3262: @*/
3263: PetscErrorCode SNESSetConvergenceHistory(SNES snes,PetscReal a[],PetscInt its[],PetscInt na,PetscBool reset)
3264: {
3271: if (na == PETSC_DECIDE || na == PETSC_DEFAULT || !a) {
3272: if (na == PETSC_DECIDE || na == PETSC_DEFAULT) na = 1000;
3273: PetscMalloc(na*sizeof(PetscReal),&a);
3274: PetscMalloc(na*sizeof(PetscInt),&its);
3275: snes->conv_malloc = PETSC_TRUE;
3276: }
3277: snes->conv_hist = a;
3278: snes->conv_hist_its = its;
3279: snes->conv_hist_max = na;
3280: snes->conv_hist_len = 0;
3281: snes->conv_hist_reset = reset;
3282: return(0);
3283: }
3285: #if defined(PETSC_HAVE_MATLAB_ENGINE)
3286: #include <engine.h> /* MATLAB include file */
3287: #include <mex.h> /* MATLAB include file */
3288: EXTERN_C_BEGIN
3291: mxArray *SNESGetConvergenceHistoryMatlab(SNES snes)
3292: {
3293: mxArray *mat;
3294: PetscInt i;
3295: PetscReal *ar;
3298: mat = mxCreateDoubleMatrix(snes->conv_hist_len,1,mxREAL);
3299: ar = (PetscReal*) mxGetData(mat);
3300: for (i=0; i<snes->conv_hist_len; i++) {
3301: ar[i] = snes->conv_hist[i];
3302: }
3303: PetscFunctionReturn(mat);
3304: }
3305: EXTERN_C_END
3306: #endif
3311: /*@C
3312: SNESGetConvergenceHistory - Gets the array used to hold the convergence history.
3314: Not Collective
3316: Input Parameter:
3317: . snes - iterative context obtained from SNESCreate()
3319: Output Parameters:
3320: . a - array to hold history
3321: . its - integer array holds the number of linear iterations (or
3322: negative if not converged) for each solve.
3323: - na - size of a and its
3325: Notes:
3326: The calling sequence for this routine in Fortran is
3327: $ call SNESGetConvergenceHistory(SNES snes, integer na, integer ierr)
3329: This routine is useful, e.g., when running a code for purposes
3330: of accurate performance monitoring, when no I/O should be done
3331: during the section of code that is being timed.
3333: Level: intermediate
3335: .keywords: SNES, get, convergence, history
3337: .seealso: SNESSetConvergencHistory()
3339: @*/
3340: PetscErrorCode SNESGetConvergenceHistory(SNES snes,PetscReal *a[],PetscInt *its[],PetscInt *na)
3341: {
3344: if (a) *a = snes->conv_hist;
3345: if (its) *its = snes->conv_hist_its;
3346: if (na) *na = snes->conv_hist_len;
3347: return(0);
3348: }
3352: /*@C
3353: SNESSetUpdate - Sets the general-purpose update function called
3354: at the beginning of every iteration of the nonlinear solve. Specifically
3355: it is called just before the Jacobian is "evaluated".
3357: Logically Collective on SNES
3359: Input Parameters:
3360: . snes - The nonlinear solver context
3361: . func - The function
3363: Calling sequence of func:
3364: . func (SNES snes, PetscInt step);
3366: . step - The current step of the iteration
3368: Level: advanced
3370: 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()
3371: This is not used by most users.
3373: .keywords: SNES, update
3375: .seealso SNESDefaultUpdate(), SNESSetJacobian(), SNESSolve()
3376: @*/
3377: PetscErrorCode SNESSetUpdate(SNES snes, PetscErrorCode (*func)(SNES, PetscInt))
3378: {
3381: snes->ops->update = func;
3382: return(0);
3383: }
3387: /*@
3388: SNESDefaultUpdate - The default update function which does nothing.
3390: Not collective
3392: Input Parameters:
3393: . snes - The nonlinear solver context
3394: . step - The current step of the iteration
3396: Level: intermediate
3398: .keywords: SNES, update
3399: .seealso SNESSetUpdate(), SNESDefaultRhsBC(), SNESDefaultShortolutionBC()
3400: @*/
3401: PetscErrorCode SNESDefaultUpdate(SNES snes, PetscInt step)
3402: {
3404: return(0);
3405: }
3409: /*
3410: SNESScaleStep_Private - Scales a step so that its length is less than the
3411: positive parameter delta.
3413: Input Parameters:
3414: + snes - the SNES context
3415: . y - approximate solution of linear system
3416: . fnorm - 2-norm of current function
3417: - delta - trust region size
3419: Output Parameters:
3420: + gpnorm - predicted function norm at the new point, assuming local
3421: linearization. The value is zero if the step lies within the trust
3422: region, and exceeds zero otherwise.
3423: - ynorm - 2-norm of the step
3425: Note:
3426: For non-trust region methods such as SNESLS, the parameter delta
3427: is set to be the maximum allowable step size.
3429: .keywords: SNES, nonlinear, scale, step
3430: */
3431: PetscErrorCode SNESScaleStep_Private(SNES snes,Vec y,PetscReal *fnorm,PetscReal *delta,PetscReal *gpnorm,PetscReal *ynorm)
3432: {
3433: PetscReal nrm;
3434: PetscScalar cnorm;
3442: VecNorm(y,NORM_2,&nrm);
3443: if (nrm > *delta) {
3444: nrm = *delta/nrm;
3445: *gpnorm = (1.0 - nrm)*(*fnorm);
3446: cnorm = nrm;
3447: VecScale(y,cnorm);
3448: *ynorm = *delta;
3449: } else {
3450: *gpnorm = 0.0;
3451: *ynorm = nrm;
3452: }
3453: return(0);
3454: }
3458: /*@C
3459: SNESSolve - Solves a nonlinear system F(x) = b.
3460: Call SNESSolve() after calling SNESCreate() and optional routines of the form SNESSetXXX().
3462: Collective on SNES
3464: Input Parameters:
3465: + snes - the SNES context
3466: . b - the constant part of the equation F(x) = b, or PETSC_NULL to use zero.
3467: - x - the solution vector.
3469: Notes:
3470: The user should initialize the vector,x, with the initial guess
3471: for the nonlinear solve prior to calling SNESSolve. In particular,
3472: to employ an initial guess of zero, the user should explicitly set
3473: this vector to zero by calling VecSet().
3475: Level: beginner
3477: .keywords: SNES, nonlinear, solve
3479: .seealso: SNESCreate(), SNESDestroy(), SNESSetFunction(), SNESSetJacobian(), SNESSetGridSequence(), SNESGetSolution()
3480: @*/
3481: PetscErrorCode SNESSolve(SNES snes,Vec b,Vec x)
3482: {
3484: PetscBool flg;
3485: char filename[PETSC_MAX_PATH_LEN];
3486: PetscViewer viewer;
3487: PetscInt grid;
3488: Vec xcreated = PETSC_NULL;
3489: DM dm;
3498: if (!x) {
3499: SNESGetDM(snes,&dm);
3500: DMCreateGlobalVector(dm,&xcreated);
3501: x = xcreated;
3502: }
3504: for (grid=0; grid<snes->gridsequence; grid++) {PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm));}
3505: for (grid=0; grid<snes->gridsequence+1; grid++) {
3507: /* set solution vector */
3508: if (!grid) {PetscObjectReference((PetscObject)x);}
3509: VecDestroy(&snes->vec_sol);
3510: snes->vec_sol = x;
3511: SNESGetDM(snes,&dm);
3513: /* set affine vector if provided */
3514: if (b) { PetscObjectReference((PetscObject)b); }
3515: VecDestroy(&snes->vec_rhs);
3516: snes->vec_rhs = b;
3518: SNESSetUp(snes);
3520: if (!grid) {
3521: if (snes->ops->computeinitialguess) {
3522: (*snes->ops->computeinitialguess)(snes,snes->vec_sol,snes->initialguessP);
3523: } else if (snes->dm) {
3524: PetscBool ig;
3525: DMHasInitialGuess(snes->dm,&ig);
3526: if (ig) {
3527: DMComputeInitialGuess(snes->dm,snes->vec_sol);
3528: }
3529: }
3530: }
3532: if (snes->conv_hist_reset) snes->conv_hist_len = 0;
3533: snes->nfuncs = 0; snes->linear_its = 0; snes->numFailures = 0;
3535: PetscLogEventBegin(SNES_Solve,snes,0,0,0);
3536: (*snes->ops->solve)(snes);
3537: PetscLogEventEnd(SNES_Solve,snes,0,0,0);
3538: if (snes->domainerror){
3539: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
3540: snes->domainerror = PETSC_FALSE;
3541: }
3542: if (!snes->reason) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
3544: flg = PETSC_FALSE;
3545: PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_test_local_min",&flg,PETSC_NULL);
3546: if (flg && !PetscPreLoadingOn) { SNESTestLocalMin(snes); }
3547: if (snes->printreason) {
3548: PetscViewerASCIIAddTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),((PetscObject)snes)->tablevel);
3549: if (snes->reason > 0) {
3550: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),"Nonlinear solve converged due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);
3551: } else {
3552: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),"Nonlinear solve did not converge due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);
3553: }
3554: PetscViewerASCIISubtractTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),((PetscObject)snes)->tablevel);
3555: }
3557: if (snes->errorifnotconverged && snes->reason < 0) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged");
3558: if (grid < snes->gridsequence) {
3559: DM fine;
3560: Vec xnew;
3561: Mat interp;
3563: DMRefine(snes->dm,((PetscObject)snes)->comm,&fine);
3564: if (!fine) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_INCOMP,"DMRefine() did not perform any refinement, cannot continue grid sequencing");
3565: DMCreateInterpolation(snes->dm,fine,&interp,PETSC_NULL);
3566: DMCreateGlobalVector(fine,&xnew);
3567: MatInterpolate(interp,x,xnew);
3568: DMInterpolate(snes->dm,interp,fine);
3569: MatDestroy(&interp);
3570: x = xnew;
3572: SNESReset(snes);
3573: SNESSetDM(snes,fine);
3574: DMDestroy(&fine);
3575: PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm));
3576: }
3577: }
3578: /* monitoring and viewing */
3579: flg = PETSC_FALSE;
3580: PetscOptionsGetString(((PetscObject)snes)->prefix,"-snes_view",filename,PETSC_MAX_PATH_LEN,&flg);
3581: if (flg && !PetscPreLoadingOn) {
3582: PetscViewerASCIIOpen(((PetscObject)snes)->comm,filename,&viewer);
3583: SNESView(snes,viewer);
3584: PetscViewerDestroy(&viewer);
3585: }
3587: flg = PETSC_FALSE;
3588: PetscOptionsGetString(((PetscObject)snes)->prefix,"-snes_view_solution_vtk",filename,PETSC_MAX_PATH_LEN,&flg);
3589: if (flg) {
3590: PetscViewer viewer;
3591: PetscViewerCreate(((PetscObject)snes)->comm,&viewer);
3592: PetscViewerSetType(viewer,PETSCVIEWERVTK);
3593: PetscViewerFileSetName(viewer,filename);
3594: VecView(snes->vec_sol,viewer);
3595: PetscViewerDestroy(&viewer);
3596: }
3598: VecDestroy(&xcreated);
3599: return(0);
3600: }
3602: /* --------- Internal routines for SNES Package --------- */
3606: /*@C
3607: SNESSetType - Sets the method for the nonlinear solver.
3609: Collective on SNES
3611: Input Parameters:
3612: + snes - the SNES context
3613: - type - a known method
3615: Options Database Key:
3616: . -snes_type <type> - Sets the method; use -help for a list
3617: of available methods (for instance, ls or tr)
3619: Notes:
3620: See "petsc/include/petscsnes.h" for available methods (for instance)
3621: + SNESLS - Newton's method with line search
3622: (systems of nonlinear equations)
3623: . SNESTR - Newton's method with trust region
3624: (systems of nonlinear equations)
3626: Normally, it is best to use the SNESSetFromOptions() command and then
3627: set the SNES solver type from the options database rather than by using
3628: this routine. Using the options database provides the user with
3629: maximum flexibility in evaluating the many nonlinear solvers.
3630: The SNESSetType() routine is provided for those situations where it
3631: is necessary to set the nonlinear solver independently of the command
3632: line or options database. This might be the case, for example, when
3633: the choice of solver changes during the execution of the program,
3634: and the user's application is taking responsibility for choosing the
3635: appropriate method.
3637: Level: intermediate
3639: .keywords: SNES, set, type
3641: .seealso: SNESType, SNESCreate()
3643: @*/
3644: PetscErrorCode SNESSetType(SNES snes,const SNESType type)
3645: {
3646: PetscErrorCode ierr,(*r)(SNES);
3647: PetscBool match;
3653: PetscObjectTypeCompare((PetscObject)snes,type,&match);
3654: if (match) return(0);
3656: PetscFListFind(SNESList,((PetscObject)snes)->comm,type,PETSC_TRUE,(void (**)(void)) &r);
3657: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested SNES type %s",type);
3658: /* Destroy the previous private SNES context */
3659: if (snes->ops->destroy) {
3660: (*(snes)->ops->destroy)(snes);
3661: snes->ops->destroy = PETSC_NULL;
3662: }
3663: /* Reinitialize function pointers in SNESOps structure */
3664: snes->ops->setup = 0;
3665: snes->ops->solve = 0;
3666: snes->ops->view = 0;
3667: snes->ops->setfromoptions = 0;
3668: snes->ops->destroy = 0;
3669: /* Call the SNESCreate_XXX routine for this particular Nonlinear solver */
3670: snes->setupcalled = PETSC_FALSE;
3671: PetscObjectChangeTypeName((PetscObject)snes,type);
3672: (*r)(snes);
3673: #if defined(PETSC_HAVE_AMS)
3674: if (PetscAMSPublishAll) {
3675: PetscObjectAMSPublish((PetscObject)snes);
3676: }
3677: #endif
3678: return(0);
3679: }
3682: /* --------------------------------------------------------------------- */
3685: /*@
3686: SNESRegisterDestroy - Frees the list of nonlinear solvers that were
3687: registered by SNESRegisterDynamic().
3689: Not Collective
3691: Level: advanced
3693: .keywords: SNES, nonlinear, register, destroy
3695: .seealso: SNESRegisterAll(), SNESRegisterAll()
3696: @*/
3697: PetscErrorCode SNESRegisterDestroy(void)
3698: {
3702: PetscFListDestroy(&SNESList);
3703: SNESRegisterAllCalled = PETSC_FALSE;
3704: return(0);
3705: }
3709: /*@C
3710: SNESGetType - Gets the SNES method type and name (as a string).
3712: Not Collective
3714: Input Parameter:
3715: . snes - nonlinear solver context
3717: Output Parameter:
3718: . type - SNES method (a character string)
3720: Level: intermediate
3722: .keywords: SNES, nonlinear, get, type, name
3723: @*/
3724: PetscErrorCode SNESGetType(SNES snes,const SNESType *type)
3725: {
3729: *type = ((PetscObject)snes)->type_name;
3730: return(0);
3731: }
3735: /*@
3736: SNESGetSolution - Returns the vector where the approximate solution is
3737: stored. This is the fine grid solution when using SNESSetGridSequence().
3739: Not Collective, but Vec is parallel if SNES is parallel
3741: Input Parameter:
3742: . snes - the SNES context
3744: Output Parameter:
3745: . x - the solution
3747: Level: intermediate
3749: .keywords: SNES, nonlinear, get, solution
3751: .seealso: SNESGetSolutionUpdate(), SNESGetFunction()
3752: @*/
3753: PetscErrorCode SNESGetSolution(SNES snes,Vec *x)
3754: {
3758: *x = snes->vec_sol;
3759: return(0);
3760: }
3764: /*@
3765: SNESGetSolutionUpdate - Returns the vector where the solution update is
3766: stored.
3768: Not Collective, but Vec is parallel if SNES is parallel
3770: Input Parameter:
3771: . snes - the SNES context
3773: Output Parameter:
3774: . x - the solution update
3776: Level: advanced
3778: .keywords: SNES, nonlinear, get, solution, update
3780: .seealso: SNESGetSolution(), SNESGetFunction()
3781: @*/
3782: PetscErrorCode SNESGetSolutionUpdate(SNES snes,Vec *x)
3783: {
3787: *x = snes->vec_sol_update;
3788: return(0);
3789: }
3793: /*@C
3794: SNESGetFunction - Returns the vector where the function is stored.
3796: Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet.
3798: Input Parameter:
3799: . snes - the SNES context
3801: Output Parameter:
3802: + r - the function (or PETSC_NULL)
3803: . func - the function (or PETSC_NULL)
3804: - ctx - the function context (or PETSC_NULL)
3806: Level: advanced
3808: .keywords: SNES, nonlinear, get, function
3810: .seealso: SNESSetFunction(), SNESGetSolution()
3811: @*/
3812: PetscErrorCode SNESGetFunction(SNES snes,Vec *r,PetscErrorCode (**func)(SNES,Vec,Vec,void*),void **ctx)
3813: {
3815: DM dm;
3819: if (r) {
3820: if (!snes->vec_func) {
3821: if (snes->vec_rhs) {
3822: VecDuplicate(snes->vec_rhs,&snes->vec_func);
3823: } else if (snes->vec_sol) {
3824: VecDuplicate(snes->vec_sol,&snes->vec_func);
3825: } else if (snes->dm) {
3826: DMCreateGlobalVector(snes->dm,&snes->vec_func);
3827: }
3828: }
3829: *r = snes->vec_func;
3830: }
3831: SNESGetDM(snes,&dm);
3832: DMSNESGetFunction(dm,func,ctx);
3833: return(0);
3834: }
3836: /*@C
3837: SNESGetGS - Returns the GS function and context.
3839: Input Parameter:
3840: . snes - the SNES context
3842: Output Parameter:
3843: + gsfunc - the function (or PETSC_NULL)
3844: - ctx - the function context (or PETSC_NULL)
3846: Level: advanced
3848: .keywords: SNES, nonlinear, get, function
3850: .seealso: SNESSetGS(), SNESGetFunction()
3851: @*/
3855: PetscErrorCode SNESGetGS (SNES snes, PetscErrorCode(**func)(SNES, Vec, Vec, void*), void ** ctx)
3856: {
3858: DM dm;
3862: SNESGetDM(snes,&dm);
3863: DMSNESGetGS(dm,func,ctx);
3864: return(0);
3865: }
3869: /*@C
3870: SNESSetOptionsPrefix - Sets the prefix used for searching for all
3871: SNES options in the database.
3873: Logically Collective on SNES
3875: Input Parameter:
3876: + snes - the SNES context
3877: - prefix - the prefix to prepend to all option names
3879: Notes:
3880: A hyphen (-) must NOT be given at the beginning of the prefix name.
3881: The first character of all runtime options is AUTOMATICALLY the hyphen.
3883: Level: advanced
3885: .keywords: SNES, set, options, prefix, database
3887: .seealso: SNESSetFromOptions()
3888: @*/
3889: PetscErrorCode SNESSetOptionsPrefix(SNES snes,const char prefix[])
3890: {
3895: PetscObjectSetOptionsPrefix((PetscObject)snes,prefix);
3896: if (!snes->ksp) {SNESGetKSP(snes,&snes->ksp);}
3897: if (snes->linesearch) {
3898: SNESGetSNESLineSearch(snes,&snes->linesearch);
3899: PetscObjectSetOptionsPrefix((PetscObject)snes->linesearch,prefix);
3900: }
3901: KSPSetOptionsPrefix(snes->ksp,prefix);
3902: return(0);
3903: }
3907: /*@C
3908: SNESAppendOptionsPrefix - Appends to the prefix used for searching for all
3909: SNES options in the database.
3911: Logically Collective on SNES
3913: Input Parameters:
3914: + snes - the SNES context
3915: - prefix - the prefix to prepend to all option names
3917: Notes:
3918: A hyphen (-) must NOT be given at the beginning of the prefix name.
3919: The first character of all runtime options is AUTOMATICALLY the hyphen.
3921: Level: advanced
3923: .keywords: SNES, append, options, prefix, database
3925: .seealso: SNESGetOptionsPrefix()
3926: @*/
3927: PetscErrorCode SNESAppendOptionsPrefix(SNES snes,const char prefix[])
3928: {
3930:
3933: PetscObjectAppendOptionsPrefix((PetscObject)snes,prefix);
3934: if (!snes->ksp) {SNESGetKSP(snes,&snes->ksp);}
3935: if (snes->linesearch) {
3936: SNESGetSNESLineSearch(snes,&snes->linesearch);
3937: PetscObjectAppendOptionsPrefix((PetscObject)snes->linesearch,prefix);
3938: }
3939: KSPAppendOptionsPrefix(snes->ksp,prefix);
3940: return(0);
3941: }
3945: /*@C
3946: SNESGetOptionsPrefix - Sets the prefix used for searching for all
3947: SNES options in the database.
3949: Not Collective
3951: Input Parameter:
3952: . snes - the SNES context
3954: Output Parameter:
3955: . prefix - pointer to the prefix string used
3957: Notes: On the fortran side, the user should pass in a string 'prefix' of
3958: sufficient length to hold the prefix.
3960: Level: advanced
3962: .keywords: SNES, get, options, prefix, database
3964: .seealso: SNESAppendOptionsPrefix()
3965: @*/
3966: PetscErrorCode SNESGetOptionsPrefix(SNES snes,const char *prefix[])
3967: {
3972: PetscObjectGetOptionsPrefix((PetscObject)snes,prefix);
3973: return(0);
3974: }
3979: /*@C
3980: SNESRegister - See SNESRegisterDynamic()
3982: Level: advanced
3983: @*/
3984: PetscErrorCode SNESRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(SNES))
3985: {
3986: char fullname[PETSC_MAX_PATH_LEN];
3990: PetscFListConcat(path,name,fullname);
3991: PetscFListAdd(&SNESList,sname,fullname,(void (*)(void))function);
3992: return(0);
3993: }
3997: PetscErrorCode SNESTestLocalMin(SNES snes)
3998: {
4000: PetscInt N,i,j;
4001: Vec u,uh,fh;
4002: PetscScalar value;
4003: PetscReal norm;
4006: SNESGetSolution(snes,&u);
4007: VecDuplicate(u,&uh);
4008: VecDuplicate(u,&fh);
4010: /* currently only works for sequential */
4011: PetscPrintf(PETSC_COMM_WORLD,"Testing FormFunction() for local min\n");
4012: VecGetSize(u,&N);
4013: for (i=0; i<N; i++) {
4014: VecCopy(u,uh);
4015: PetscPrintf(PETSC_COMM_WORLD,"i = %D\n",i);
4016: for (j=-10; j<11; j++) {
4017: value = PetscSign(j)*exp(PetscAbs(j)-10.0);
4018: VecSetValue(uh,i,value,ADD_VALUES);
4019: SNESComputeFunction(snes,uh,fh);
4020: VecNorm(fh,NORM_2,&norm);
4021: PetscPrintf(PETSC_COMM_WORLD," j norm %D %18.16e\n",j,norm);
4022: value = -value;
4023: VecSetValue(uh,i,value,ADD_VALUES);
4024: }
4025: }
4026: VecDestroy(&uh);
4027: VecDestroy(&fh);
4028: return(0);
4029: }
4033: /*@
4034: SNESKSPSetUseEW - Sets SNES use Eisenstat-Walker method for
4035: computing relative tolerance for linear solvers within an inexact
4036: Newton method.
4038: Logically Collective on SNES
4040: Input Parameters:
4041: + snes - SNES context
4042: - flag - PETSC_TRUE or PETSC_FALSE
4044: Options Database:
4045: + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence
4046: . -snes_ksp_ew_version ver - version of Eisenstat-Walker method
4047: . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0
4048: . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax
4049: . -snes_ksp_ew_gamma <gamma> - Sets gamma
4050: . -snes_ksp_ew_alpha <alpha> - Sets alpha
4051: . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2
4052: - -snes_ksp_ew_threshold <threshold> - Sets threshold
4054: Notes:
4055: Currently, the default is to use a constant relative tolerance for
4056: the inner linear solvers. Alternatively, one can use the
4057: Eisenstat-Walker method, where the relative convergence tolerance
4058: is reset at each Newton iteration according progress of the nonlinear
4059: solver.
4061: Level: advanced
4063: Reference:
4064: S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an
4065: inexact Newton method", SISC 17 (1), pp.16-32, 1996.
4067: .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton
4069: .seealso: SNESKSPGetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW()
4070: @*/
4071: PetscErrorCode SNESKSPSetUseEW(SNES snes,PetscBool flag)
4072: {
4076: snes->ksp_ewconv = flag;
4077: return(0);
4078: }
4082: /*@
4083: SNESKSPGetUseEW - Gets if SNES is using Eisenstat-Walker method
4084: for computing relative tolerance for linear solvers within an
4085: inexact Newton method.
4087: Not Collective
4089: Input Parameter:
4090: . snes - SNES context
4092: Output Parameter:
4093: . flag - PETSC_TRUE or PETSC_FALSE
4095: Notes:
4096: Currently, the default is to use a constant relative tolerance for
4097: the inner linear solvers. Alternatively, one can use the
4098: Eisenstat-Walker method, where the relative convergence tolerance
4099: is reset at each Newton iteration according progress of the nonlinear
4100: solver.
4102: Level: advanced
4104: Reference:
4105: S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an
4106: inexact Newton method", SISC 17 (1), pp.16-32, 1996.
4108: .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton
4110: .seealso: SNESKSPSetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW()
4111: @*/
4112: PetscErrorCode SNESKSPGetUseEW(SNES snes, PetscBool *flag)
4113: {
4117: *flag = snes->ksp_ewconv;
4118: return(0);
4119: }
4123: /*@
4124: SNESKSPSetParametersEW - Sets parameters for Eisenstat-Walker
4125: convergence criteria for the linear solvers within an inexact
4126: Newton method.
4128: Logically Collective on SNES
4129:
4130: Input Parameters:
4131: + snes - SNES context
4132: . version - version 1, 2 (default is 2) or 3
4133: . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1)
4134: . rtol_max - maximum relative tolerance (0 <= rtol_max < 1)
4135: . gamma - multiplicative factor for version 2 rtol computation
4136: (0 <= gamma2 <= 1)
4137: . alpha - power for version 2 rtol computation (1 < alpha <= 2)
4138: . alpha2 - power for safeguard
4139: - threshold - threshold for imposing safeguard (0 < threshold < 1)
4141: Note:
4142: Version 3 was contributed by Luis Chacon, June 2006.
4144: Use PETSC_DEFAULT to retain the default for any of the parameters.
4146: Level: advanced
4148: Reference:
4149: S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an
4150: inexact Newton method", Utah State University Math. Stat. Dept. Res.
4151: Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput.
4153: .keywords: SNES, KSP, Eisenstat, Walker, set, parameters
4155: .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPGetParametersEW()
4156: @*/
4157: PetscErrorCode SNESKSPSetParametersEW(SNES snes,PetscInt version,PetscReal rtol_0,PetscReal rtol_max,
4158: PetscReal gamma,PetscReal alpha,PetscReal alpha2,PetscReal threshold)
4159: {
4160: SNESKSPEW *kctx;
4163: kctx = (SNESKSPEW*)snes->kspconvctx;
4164: if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing");
4173: if (version != PETSC_DEFAULT) kctx->version = version;
4174: if (rtol_0 != PETSC_DEFAULT) kctx->rtol_0 = rtol_0;
4175: if (rtol_max != PETSC_DEFAULT) kctx->rtol_max = rtol_max;
4176: if (gamma != PETSC_DEFAULT) kctx->gamma = gamma;
4177: if (alpha != PETSC_DEFAULT) kctx->alpha = alpha;
4178: if (alpha2 != PETSC_DEFAULT) kctx->alpha2 = alpha2;
4179: if (threshold != PETSC_DEFAULT) kctx->threshold = threshold;
4180:
4181: if (kctx->version < 1 || kctx->version > 3) {
4182: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 and 3 are supported: %D",kctx->version);
4183: }
4184: if (kctx->rtol_0 < 0.0 || kctx->rtol_0 >= 1.0) {
4185: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= rtol_0 < 1.0: %G",kctx->rtol_0);
4186: }
4187: if (kctx->rtol_max < 0.0 || kctx->rtol_max >= 1.0) {
4188: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= rtol_max (%G) < 1.0\n",kctx->rtol_max);
4189: }
4190: if (kctx->gamma < 0.0 || kctx->gamma > 1.0) {
4191: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= gamma (%G) <= 1.0\n",kctx->gamma);
4192: }
4193: if (kctx->alpha <= 1.0 || kctx->alpha > 2.0) {
4194: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"1.0 < alpha (%G) <= 2.0\n",kctx->alpha);
4195: }
4196: if (kctx->threshold <= 0.0 || kctx->threshold >= 1.0) {
4197: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 < threshold (%G) < 1.0\n",kctx->threshold);
4198: }
4199: return(0);
4200: }
4204: /*@
4205: SNESKSPGetParametersEW - Gets parameters for Eisenstat-Walker
4206: convergence criteria for the linear solvers within an inexact
4207: Newton method.
4209: Not Collective
4210:
4211: Input Parameters:
4212: snes - SNES context
4214: Output Parameters:
4215: + version - version 1, 2 (default is 2) or 3
4216: . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1)
4217: . rtol_max - maximum relative tolerance (0 <= rtol_max < 1)
4218: . gamma - multiplicative factor for version 2 rtol computation
4219: (0 <= gamma2 <= 1)
4220: . alpha - power for version 2 rtol computation (1 < alpha <= 2)
4221: . alpha2 - power for safeguard
4222: - threshold - threshold for imposing safeguard (0 < threshold < 1)
4224: Level: advanced
4226: .keywords: SNES, KSP, Eisenstat, Walker, get, parameters
4228: .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPSetParametersEW()
4229: @*/
4230: PetscErrorCode SNESKSPGetParametersEW(SNES snes,PetscInt *version,PetscReal *rtol_0,PetscReal *rtol_max,
4231: PetscReal *gamma,PetscReal *alpha,PetscReal *alpha2,PetscReal *threshold)
4232: {
4233: SNESKSPEW *kctx;
4236: kctx = (SNESKSPEW*)snes->kspconvctx;
4237: if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing");
4238: if(version) *version = kctx->version;
4239: if(rtol_0) *rtol_0 = kctx->rtol_0;
4240: if(rtol_max) *rtol_max = kctx->rtol_max;
4241: if(gamma) *gamma = kctx->gamma;
4242: if(alpha) *alpha = kctx->alpha;
4243: if(alpha2) *alpha2 = kctx->alpha2;
4244: if(threshold) *threshold = kctx->threshold;
4245: return(0);
4246: }
4250: static PetscErrorCode SNESKSPEW_PreSolve(SNES snes, KSP ksp, Vec b, Vec x)
4251: {
4253: SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx;
4254: PetscReal rtol=PETSC_DEFAULT,stol;
4257: if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context exists");
4258: if (!snes->iter) { /* first time in, so use the original user rtol */
4259: rtol = kctx->rtol_0;
4260: } else {
4261: if (kctx->version == 1) {
4262: rtol = (snes->norm - kctx->lresid_last)/kctx->norm_last;
4263: if (rtol < 0.0) rtol = -rtol;
4264: stol = pow(kctx->rtol_last,kctx->alpha2);
4265: if (stol > kctx->threshold) rtol = PetscMax(rtol,stol);
4266: } else if (kctx->version == 2) {
4267: rtol = kctx->gamma * pow(snes->norm/kctx->norm_last,kctx->alpha);
4268: stol = kctx->gamma * pow(kctx->rtol_last,kctx->alpha);
4269: if (stol > kctx->threshold) rtol = PetscMax(rtol,stol);
4270: } else if (kctx->version == 3) {/* contributed by Luis Chacon, June 2006. */
4271: rtol = kctx->gamma * pow(snes->norm/kctx->norm_last,kctx->alpha);
4272: /* safeguard: avoid sharp decrease of rtol */
4273: stol = kctx->gamma*pow(kctx->rtol_last,kctx->alpha);
4274: stol = PetscMax(rtol,stol);
4275: rtol = PetscMin(kctx->rtol_0,stol);
4276: /* safeguard: avoid oversolving */
4277: stol = kctx->gamma*(snes->ttol)/snes->norm;
4278: stol = PetscMax(rtol,stol);
4279: rtol = PetscMin(kctx->rtol_0,stol);
4280: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 or 3 are supported: %D",kctx->version);
4281: }
4282: /* safeguard: avoid rtol greater than one */
4283: rtol = PetscMin(rtol,kctx->rtol_max);
4284: KSPSetTolerances(ksp,rtol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
4285: PetscInfo3(snes,"iter %D, Eisenstat-Walker (version %D) KSP rtol=%G\n",snes->iter,kctx->version,rtol);
4286: return(0);
4287: }
4291: static PetscErrorCode SNESKSPEW_PostSolve(SNES snes, KSP ksp, Vec b, Vec x)
4292: {
4294: SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx;
4295: PCSide pcside;
4296: Vec lres;
4299: if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context exists");
4300: KSPGetTolerances(ksp,&kctx->rtol_last,0,0,0);
4301: SNESGetFunctionNorm(snes,&kctx->norm_last);
4302: if (kctx->version == 1) {
4303: KSPGetPCSide(ksp,&pcside);
4304: if (pcside == PC_RIGHT) { /* XXX Should we also test KSP_UNPRECONDITIONED_NORM ? */
4305: /* KSP residual is true linear residual */
4306: KSPGetResidualNorm(ksp,&kctx->lresid_last);
4307: } else {
4308: /* KSP residual is preconditioned residual */
4309: /* compute true linear residual norm */
4310: VecDuplicate(b,&lres);
4311: MatMult(snes->jacobian,x,lres);
4312: VecAYPX(lres,-1.0,b);
4313: VecNorm(lres,NORM_2,&kctx->lresid_last);
4314: VecDestroy(&lres);
4315: }
4316: }
4317: return(0);
4318: }
4322: PetscErrorCode SNES_KSPSolve(SNES snes, KSP ksp, Vec b, Vec x)
4323: {
4327: if (snes->ksp_ewconv) { SNESKSPEW_PreSolve(snes,ksp,b,x); }
4328: KSPSolve(ksp,b,x);
4329: if (snes->ksp_ewconv) { SNESKSPEW_PostSolve(snes,ksp,b,x); }
4330: return(0);
4331: }
4335: /*@
4336: SNESSetDM - Sets the DM that may be used by some preconditioners
4338: Logically Collective on SNES
4340: Input Parameters:
4341: + snes - the preconditioner context
4342: - dm - the dm
4344: Level: intermediate
4347: .seealso: SNESGetDM(), KSPSetDM(), KSPGetDM()
4348: @*/
4349: PetscErrorCode SNESSetDM(SNES snes,DM dm)
4350: {
4352: KSP ksp;
4353: SNESDM sdm;
4357: if (dm) {PetscObjectReference((PetscObject)dm);}
4358: if (snes->dm) { /* Move the SNESDM context over to the new DM unless the new DM already has one */
4359: PetscContainer oldcontainer,container;
4360: PetscObjectQuery((PetscObject)snes->dm,"SNESDM",(PetscObject*)&oldcontainer);
4361: PetscObjectQuery((PetscObject)dm,"SNESDM",(PetscObject*)&container);
4362: if (oldcontainer && !container) {
4363: DMSNESCopyContext(snes->dm,dm);
4364: DMSNESGetContext(snes->dm,&sdm);
4365: if (sdm->originaldm == snes->dm) { /* Grant write privileges to the replacement DM */
4366: sdm->originaldm = dm;
4367: }
4368: }
4369: DMDestroy(&snes->dm);
4370: }
4371: snes->dm = dm;
4372: SNESGetKSP(snes,&ksp);
4373: KSPSetDM(ksp,dm);
4374: KSPSetDMActive(ksp,PETSC_FALSE);
4375: if (snes->pc) {
4376: SNESSetDM(snes->pc, snes->dm);
4377: }
4378: return(0);
4379: }
4383: /*@
4384: SNESGetDM - Gets the DM that may be used by some preconditioners
4386: Not Collective but DM obtained is parallel on SNES
4388: Input Parameter:
4389: . snes - the preconditioner context
4391: Output Parameter:
4392: . dm - the dm
4394: Level: intermediate
4397: .seealso: SNESSetDM(), KSPSetDM(), KSPGetDM()
4398: @*/
4399: PetscErrorCode SNESGetDM(SNES snes,DM *dm)
4400: {
4405: if (!snes->dm) {
4406: DMShellCreate(((PetscObject)snes)->comm,&snes->dm);
4407: }
4408: *dm = snes->dm;
4409: return(0);
4410: }
4414: /*@
4415: SNESSetPC - Sets the nonlinear preconditioner to be used.
4417: Collective on SNES
4419: Input Parameters:
4420: + snes - iterative context obtained from SNESCreate()
4421: - pc - the preconditioner object
4423: Notes:
4424: Use SNESGetPC() to retrieve the preconditioner context (for example,
4425: to configure it using the API).
4427: Level: developer
4429: .keywords: SNES, set, precondition
4430: .seealso: SNESGetPC()
4431: @*/
4432: PetscErrorCode SNESSetPC(SNES snes, SNES pc)
4433: {
4440: PetscObjectReference((PetscObject) pc);
4441: SNESDestroy(&snes->pc);
4442: snes->pc = pc;
4443: PetscLogObjectParent(snes, snes->pc);
4444: return(0);
4445: }
4449: /*@
4450: SNESGetPC - Returns a pointer to the nonlinear preconditioning context set with SNESSetPC().
4452: Not Collective
4454: Input Parameter:
4455: . snes - iterative context obtained from SNESCreate()
4457: Output Parameter:
4458: . pc - preconditioner context
4460: Level: developer
4462: .keywords: SNES, get, preconditioner
4463: .seealso: SNESSetPC()
4464: @*/
4465: PetscErrorCode SNESGetPC(SNES snes, SNES *pc)
4466: {
4467: PetscErrorCode ierr;
4468: const char *optionsprefix;
4473: if (!snes->pc) {
4474: SNESCreate(((PetscObject) snes)->comm,&snes->pc);
4475: PetscObjectIncrementTabLevel((PetscObject)snes->pc,(PetscObject)snes,1);
4476: PetscLogObjectParent(snes,snes->pc);
4477: SNESGetOptionsPrefix(snes,&optionsprefix);
4478: SNESSetOptionsPrefix(snes->pc,optionsprefix);
4479: SNESAppendOptionsPrefix(snes->pc,"npc_");
4480: }
4481: *pc = snes->pc;
4482: return(0);
4483: }
4487: /*@
4488: SNESSetSNESLineSearch - Sets the linesearch on the SNES instance.
4490: Collective on SNES
4492: Input Parameters:
4493: + snes - iterative context obtained from SNESCreate()
4494: - linesearch - the linesearch object
4496: Notes:
4497: Use SNESGetSNESLineSearch() to retrieve the preconditioner context (for example,
4498: to configure it using the API).
4500: Level: developer
4502: .keywords: SNES, set, linesearch
4503: .seealso: SNESGetSNESLineSearch()
4504: @*/
4505: PetscErrorCode SNESSetSNESLineSearch(SNES snes, SNESLineSearch linesearch)
4506: {
4513: PetscObjectReference((PetscObject) linesearch);
4514: SNESLineSearchDestroy(&snes->linesearch);
4515: snes->linesearch = linesearch;
4516: PetscLogObjectParent(snes, snes->linesearch);
4517: return(0);
4518: }
4522: /*@C
4523: SNESGetSNESLineSearch - Returns a pointer to the line search context set with SNESSetLineSearch()
4524: or creates a default line search instance associated with the SNES and returns it.
4526: Not Collective
4528: Input Parameter:
4529: . snes - iterative context obtained from SNESCreate()
4531: Output Parameter:
4532: . linesearch - linesearch context
4534: Level: developer
4536: .keywords: SNES, get, linesearch
4537: .seealso: SNESSetSNESLineSearch()
4538: @*/
4539: PetscErrorCode SNESGetSNESLineSearch(SNES snes, SNESLineSearch *linesearch)
4540: {
4542: const char *optionsprefix;
4547: if (!snes->linesearch) {
4548: SNESGetOptionsPrefix(snes, &optionsprefix);
4549: SNESLineSearchCreate(((PetscObject) snes)->comm, &snes->linesearch);
4550: SNESLineSearchSetSNES(snes->linesearch, snes);
4551: SNESLineSearchAppendOptionsPrefix(snes->linesearch, optionsprefix);
4552: PetscObjectIncrementTabLevel((PetscObject) snes->linesearch, (PetscObject) snes, 1);
4553: PetscLogObjectParent(snes, snes->linesearch);
4554: }
4555: *linesearch = snes->linesearch;
4556: return(0);
4557: }
4559: #if defined(PETSC_HAVE_MATLAB_ENGINE)
4560: #include <mex.h>
4562: typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext;
4566: /*
4567: SNESComputeFunction_Matlab - Calls the function that has been set with
4568: SNESSetFunctionMatlab().
4570: Collective on SNES
4572: Input Parameters:
4573: + snes - the SNES context
4574: - x - input vector
4576: Output Parameter:
4577: . y - function vector, as set by SNESSetFunction()
4579: Notes:
4580: SNESComputeFunction() is typically used within nonlinear solvers
4581: implementations, so most users would not generally call this routine
4582: themselves.
4584: Level: developer
4586: .keywords: SNES, nonlinear, compute, function
4588: .seealso: SNESSetFunction(), SNESGetFunction()
4589: */
4590: PetscErrorCode SNESComputeFunction_Matlab(SNES snes,Vec x,Vec y, void *ctx)
4591: {
4592: PetscErrorCode ierr;
4593: SNESMatlabContext *sctx = (SNESMatlabContext *)ctx;
4594: int nlhs = 1,nrhs = 5;
4595: mxArray *plhs[1],*prhs[5];
4596: long long int lx = 0,ly = 0,ls = 0;
4605: /* call Matlab function in ctx with arguments x and y */
4607: PetscMemcpy(&ls,&snes,sizeof(snes));
4608: PetscMemcpy(&lx,&x,sizeof(x));
4609: PetscMemcpy(&ly,&y,sizeof(x));
4610: prhs[0] = mxCreateDoubleScalar((double)ls);
4611: prhs[1] = mxCreateDoubleScalar((double)lx);
4612: prhs[2] = mxCreateDoubleScalar((double)ly);
4613: prhs[3] = mxCreateString(sctx->funcname);
4614: prhs[4] = sctx->ctx;
4615: mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeFunctionInternal");
4616: mxGetScalar(plhs[0]);
4617: mxDestroyArray(prhs[0]);
4618: mxDestroyArray(prhs[1]);
4619: mxDestroyArray(prhs[2]);
4620: mxDestroyArray(prhs[3]);
4621: mxDestroyArray(plhs[0]);
4622: return(0);
4623: }
4628: /*
4629: SNESSetFunctionMatlab - Sets the function evaluation routine and function
4630: vector for use by the SNES routines in solving systems of nonlinear
4631: equations from MATLAB. Here the function is a string containing the name of a MATLAB function
4633: Logically Collective on SNES
4635: Input Parameters:
4636: + snes - the SNES context
4637: . r - vector to store function value
4638: - func - function evaluation routine
4640: Calling sequence of func:
4641: $ func (SNES snes,Vec x,Vec f,void *ctx);
4644: Notes:
4645: The Newton-like methods typically solve linear systems of the form
4646: $ f'(x) x = -f(x),
4647: where f'(x) denotes the Jacobian matrix and f(x) is the function.
4649: Level: beginner
4651: .keywords: SNES, nonlinear, set, function
4653: .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction()
4654: */
4655: PetscErrorCode SNESSetFunctionMatlab(SNES snes,Vec r,const char *func,mxArray *ctx)
4656: {
4657: PetscErrorCode ierr;
4658: SNESMatlabContext *sctx;
4661: /* currently sctx is memory bleed */
4662: PetscMalloc(sizeof(SNESMatlabContext),&sctx);
4663: PetscStrallocpy(func,&sctx->funcname);
4664: /*
4665: This should work, but it doesn't
4666: sctx->ctx = ctx;
4667: mexMakeArrayPersistent(sctx->ctx);
4668: */
4669: sctx->ctx = mxDuplicateArray(ctx);
4670: SNESSetFunction(snes,r,SNESComputeFunction_Matlab,sctx);
4671: return(0);
4672: }
4676: /*
4677: SNESComputeJacobian_Matlab - Calls the function that has been set with
4678: SNESSetJacobianMatlab().
4680: Collective on SNES
4682: Input Parameters:
4683: + snes - the SNES context
4684: . x - input vector
4685: . A, B - the matrices
4686: - ctx - user context
4688: Output Parameter:
4689: . flag - structure of the matrix
4691: Level: developer
4693: .keywords: SNES, nonlinear, compute, function
4695: .seealso: SNESSetFunction(), SNESGetFunction()
4696: @*/
4697: PetscErrorCode SNESComputeJacobian_Matlab(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *flag, void *ctx)
4698: {
4699: PetscErrorCode ierr;
4700: SNESMatlabContext *sctx = (SNESMatlabContext *)ctx;
4701: int nlhs = 2,nrhs = 6;
4702: mxArray *plhs[2],*prhs[6];
4703: long long int lx = 0,lA = 0,ls = 0, lB = 0;
4704:
4709: /* call Matlab function in ctx with arguments x and y */
4711: PetscMemcpy(&ls,&snes,sizeof(snes));
4712: PetscMemcpy(&lx,&x,sizeof(x));
4713: PetscMemcpy(&lA,A,sizeof(x));
4714: PetscMemcpy(&lB,B,sizeof(x));
4715: prhs[0] = mxCreateDoubleScalar((double)ls);
4716: prhs[1] = mxCreateDoubleScalar((double)lx);
4717: prhs[2] = mxCreateDoubleScalar((double)lA);
4718: prhs[3] = mxCreateDoubleScalar((double)lB);
4719: prhs[4] = mxCreateString(sctx->funcname);
4720: prhs[5] = sctx->ctx;
4721: mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeJacobianInternal");
4722: mxGetScalar(plhs[0]);
4723: *flag = (MatStructure) mxGetScalar(plhs[1]);
4724: mxDestroyArray(prhs[0]);
4725: mxDestroyArray(prhs[1]);
4726: mxDestroyArray(prhs[2]);
4727: mxDestroyArray(prhs[3]);
4728: mxDestroyArray(prhs[4]);
4729: mxDestroyArray(plhs[0]);
4730: mxDestroyArray(plhs[1]);
4731: return(0);
4732: }
4737: /*
4738: SNESSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices
4739: vector for use by the SNES routines in solving systems of nonlinear
4740: equations from MATLAB. Here the function is a string containing the name of a MATLAB function
4742: Logically Collective on SNES
4744: Input Parameters:
4745: + snes - the SNES context
4746: . A,B - Jacobian matrices
4747: . func - function evaluation routine
4748: - ctx - user context
4750: Calling sequence of func:
4751: $ flag = func (SNES snes,Vec x,Mat A,Mat B,void *ctx);
4754: Level: developer
4756: .keywords: SNES, nonlinear, set, function
4758: .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction()
4759: */
4760: PetscErrorCode SNESSetJacobianMatlab(SNES snes,Mat A,Mat B,const char *func,mxArray *ctx)
4761: {
4762: PetscErrorCode ierr;
4763: SNESMatlabContext *sctx;
4766: /* currently sctx is memory bleed */
4767: PetscMalloc(sizeof(SNESMatlabContext),&sctx);
4768: PetscStrallocpy(func,&sctx->funcname);
4769: /*
4770: This should work, but it doesn't
4771: sctx->ctx = ctx;
4772: mexMakeArrayPersistent(sctx->ctx);
4773: */
4774: sctx->ctx = mxDuplicateArray(ctx);
4775: SNESSetJacobian(snes,A,B,SNESComputeJacobian_Matlab,sctx);
4776: return(0);
4777: }
4781: /*
4782: SNESMonitor_Matlab - Calls the function that has been set with SNESMonitorSetMatlab().
4784: Collective on SNES
4786: .seealso: SNESSetFunction(), SNESGetFunction()
4787: @*/
4788: PetscErrorCode SNESMonitor_Matlab(SNES snes,PetscInt it, PetscReal fnorm, void *ctx)
4789: {
4790: PetscErrorCode ierr;
4791: SNESMatlabContext *sctx = (SNESMatlabContext *)ctx;
4792: int nlhs = 1,nrhs = 6;
4793: mxArray *plhs[1],*prhs[6];
4794: long long int lx = 0,ls = 0;
4795: Vec x=snes->vec_sol;
4796:
4800: PetscMemcpy(&ls,&snes,sizeof(snes));
4801: PetscMemcpy(&lx,&x,sizeof(x));
4802: prhs[0] = mxCreateDoubleScalar((double)ls);
4803: prhs[1] = mxCreateDoubleScalar((double)it);
4804: prhs[2] = mxCreateDoubleScalar((double)fnorm);
4805: prhs[3] = mxCreateDoubleScalar((double)lx);
4806: prhs[4] = mxCreateString(sctx->funcname);
4807: prhs[5] = sctx->ctx;
4808: mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESMonitorInternal");
4809: mxGetScalar(plhs[0]);
4810: mxDestroyArray(prhs[0]);
4811: mxDestroyArray(prhs[1]);
4812: mxDestroyArray(prhs[2]);
4813: mxDestroyArray(prhs[3]);
4814: mxDestroyArray(prhs[4]);
4815: mxDestroyArray(plhs[0]);
4816: return(0);
4817: }
4822: /*
4823: SNESMonitorSetMatlab - Sets the monitor function from MATLAB
4825: Level: developer
4827: .keywords: SNES, nonlinear, set, function
4829: .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction()
4830: */
4831: PetscErrorCode SNESMonitorSetMatlab(SNES snes,const char *func,mxArray *ctx)
4832: {
4833: PetscErrorCode ierr;
4834: SNESMatlabContext *sctx;
4837: /* currently sctx is memory bleed */
4838: PetscMalloc(sizeof(SNESMatlabContext),&sctx);
4839: PetscStrallocpy(func,&sctx->funcname);
4840: /*
4841: This should work, but it doesn't
4842: sctx->ctx = ctx;
4843: mexMakeArrayPersistent(sctx->ctx);
4844: */
4845: sctx->ctx = mxDuplicateArray(ctx);
4846: SNESMonitorSet(snes,SNESMonitor_Matlab,sctx,PETSC_NULL);
4847: return(0);
4848: }
4850: #endif