Actual source code: ls.c
2: #include <../src/snes/impls/ls/lsimpl.h>
4: /*
5: Checks if J^T F = 0 which implies we've found a local minimum of the norm of the function,
6: || F(u) ||_2 but not a zero, F(u) = 0. In the case when one cannot compute J^T F we use the fact that
7: 0 = (J^T F)^T W = F^T J W iff W not in the null space of J. Thanks for Jorge More
8: for this trick. One assumes that the probability that W is in the null space of J is very, very small.
9: */
10: static PetscErrorCode SNESNEWTONLSCheckLocalMin_Private(SNES snes, Mat A, Vec F, PetscReal fnorm, PetscBool *ismin)
11: {
12: PetscReal a1;
13: PetscBool hastranspose;
14: Vec W;
15: PetscErrorCode (*objective)(SNES, Vec, PetscReal *, void *);
17: PetscFunctionBegin;
18: *ismin = PETSC_FALSE;
19: PetscCall(SNESGetObjective(snes, &objective, NULL));
20: if (!objective) {
21: PetscCall(MatHasOperation(A, MATOP_MULT_TRANSPOSE, &hastranspose));
22: PetscCall(VecDuplicate(F, &W));
23: if (hastranspose) {
24: /* Compute || J^T F|| */
25: PetscCall(MatMultTranspose(A, F, W));
26: PetscCall(VecNorm(W, NORM_2, &a1));
27: PetscCall(PetscInfo(snes, "|| J^T F|| %14.12e near zero implies found a local minimum\n", (double)(a1 / fnorm)));
28: if (a1 / fnorm < 1.e-4) *ismin = PETSC_TRUE;
29: } else {
30: Vec work;
31: PetscScalar result;
32: PetscReal wnorm;
34: PetscCall(VecSetRandom(W, NULL));
35: PetscCall(VecNorm(W, NORM_2, &wnorm));
36: PetscCall(VecDuplicate(W, &work));
37: PetscCall(MatMult(A, W, work));
38: PetscCall(VecDot(F, work, &result));
39: PetscCall(VecDestroy(&work));
40: a1 = PetscAbsScalar(result) / (fnorm * wnorm);
41: PetscCall(PetscInfo(snes, "(F^T J random)/(|| F ||*||random|| %14.12e near zero implies found a local minimum\n", (double)a1));
42: if (a1 < 1.e-4) *ismin = PETSC_TRUE;
43: }
44: PetscCall(VecDestroy(&W));
45: }
46: PetscFunctionReturn(PETSC_SUCCESS);
47: }
49: /*
50: Checks if J^T(F - J*X) = 0
51: */
52: static PetscErrorCode SNESNEWTONLSCheckResidual_Private(SNES snes, Mat A, Vec F, Vec X)
53: {
54: PetscReal a1, a2;
55: PetscBool hastranspose;
56: PetscErrorCode (*objective)(SNES, Vec, PetscReal *, void *);
58: PetscFunctionBegin;
59: PetscCall(MatHasOperation(A, MATOP_MULT_TRANSPOSE, &hastranspose));
60: PetscCall(SNESGetObjective(snes, &objective, NULL));
61: if (hastranspose && !objective) {
62: Vec W1, W2;
64: PetscCall(VecDuplicate(F, &W1));
65: PetscCall(VecDuplicate(F, &W2));
66: PetscCall(MatMult(A, X, W1));
67: PetscCall(VecAXPY(W1, -1.0, F));
69: /* Compute || J^T W|| */
70: PetscCall(MatMultTranspose(A, W1, W2));
71: PetscCall(VecNorm(W1, NORM_2, &a1));
72: PetscCall(VecNorm(W2, NORM_2, &a2));
73: if (a1 != 0.0) PetscCall(PetscInfo(snes, "||J^T(F-Ax)||/||F-AX|| %14.12e near zero implies inconsistent rhs\n", (double)(a2 / a1)));
74: PetscCall(VecDestroy(&W1));
75: PetscCall(VecDestroy(&W2));
76: }
77: PetscFunctionReturn(PETSC_SUCCESS);
78: }
80: /*
81: This file implements a truncated Newton method with a line search,
82: for solving a system of nonlinear equations, using the KSP, Vec,
83: and Mat interfaces for linear solvers, vectors, and matrices,
84: respectively.
86: The following basic routines are required for each nonlinear solver:
87: SNESCreate_XXX() - Creates a nonlinear solver context
88: SNESSetFromOptions_XXX() - Sets runtime options
89: SNESSolve_XXX() - Solves the nonlinear system
90: SNESDestroy_XXX() - Destroys the nonlinear solver context
91: The suffix "_XXX" denotes a particular implementation, in this case
92: we use _NEWTONLS (e.g., SNESCreate_NEWTONLS, SNESSolve_NEWTONLS) for solving
93: systems of nonlinear equations with a line search (LS) method.
94: These routines are actually called via the common user interface
95: routines SNESCreate(), SNESSetFromOptions(), SNESSolve(), and
96: SNESDestroy(), so the application code interface remains identical
97: for all nonlinear solvers.
99: Another key routine is:
100: SNESSetUp_XXX() - Prepares for the use of a nonlinear solver
101: by setting data structures and options. The interface routine SNESSetUp()
102: is not usually called directly by the user, but instead is called by
103: SNESSolve() if necessary.
105: Additional basic routines are:
106: SNESView_XXX() - Prints details of runtime options that
107: have actually been used.
108: These are called by application codes via the interface routines
109: SNESView().
111: The various types of solvers (preconditioners, Krylov subspace methods,
112: nonlinear solvers, timesteppers) are all organized similarly, so the
113: above description applies to these categories also.
115: */
116: /*
117: SNESSolve_NEWTONLS - Solves a nonlinear system with a truncated Newton
118: method with a line search.
120: Input Parameter:
121: . snes - the SNES context
123: Output Parameter:
124: . outits - number of iterations until termination
126: Application Interface Routine: SNESSolve()
128: */
129: PetscErrorCode SNESSolve_NEWTONLS(SNES snes)
130: {
131: PetscInt maxits, i, lits;
132: SNESLineSearchReason lssucceed;
133: PetscReal fnorm, xnorm, ynorm;
134: Vec Y, X, F;
135: SNESLineSearch linesearch;
136: SNESConvergedReason reason;
137: #if defined(PETSC_USE_INFO)
138: PetscReal gnorm;
139: #endif
141: PetscFunctionBegin;
142: PetscCheck(!snes->xl && !snes->xu && !snes->ops->computevariablebounds, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);
144: snes->numFailures = 0;
145: snes->numLinearSolveFailures = 0;
146: snes->reason = SNES_CONVERGED_ITERATING;
148: maxits = snes->max_its; /* maximum number of iterations */
149: X = snes->vec_sol; /* solution vector */
150: F = snes->vec_func; /* residual vector */
151: Y = snes->vec_sol_update; /* newton step */
153: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
154: snes->iter = 0;
155: snes->norm = 0.0;
156: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
157: PetscCall(SNESGetLineSearch(snes, &linesearch));
159: /* compute the preconditioned function first in the case of left preconditioning with preconditioned function */
160: if (snes->npc && snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_PRECONDITIONED) {
161: PetscCall(SNESApplyNPC(snes, X, NULL, F));
162: PetscCall(SNESGetConvergedReason(snes->npc, &reason));
163: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT && reason != SNES_DIVERGED_TR_DELTA) {
164: snes->reason = SNES_DIVERGED_INNER;
165: PetscFunctionReturn(PETSC_SUCCESS);
166: }
168: PetscCall(VecNormBegin(F, NORM_2, &fnorm));
169: PetscCall(VecNormEnd(F, NORM_2, &fnorm));
170: } else {
171: if (!snes->vec_func_init_set) {
172: PetscCall(SNESComputeFunction(snes, X, F));
173: } else snes->vec_func_init_set = PETSC_FALSE;
174: }
176: PetscCall(VecNorm(F, NORM_2, &fnorm)); /* fnorm <- ||F|| */
177: SNESCheckFunctionNorm(snes, fnorm);
178: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
179: snes->norm = fnorm;
180: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
181: PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0));
182: PetscCall(SNESMonitor(snes, 0, fnorm));
184: /* test convergence */
185: PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
186: if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
188: for (i = 0; i < maxits; i++) {
189: /* Call general purpose update function */
190: PetscTryTypeMethod(snes, update, snes->iter);
192: /* apply the nonlinear preconditioner */
193: if (snes->npc) {
194: if (snes->npcside == PC_RIGHT) {
195: PetscCall(SNESSetInitialFunction(snes->npc, F));
196: PetscCall(PetscLogEventBegin(SNES_NPCSolve, snes->npc, X, snes->vec_rhs, 0));
197: PetscCall(SNESSolve(snes->npc, snes->vec_rhs, X));
198: PetscCall(PetscLogEventEnd(SNES_NPCSolve, snes->npc, X, snes->vec_rhs, 0));
199: PetscCall(SNESGetConvergedReason(snes->npc, &reason));
200: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT && reason != SNES_DIVERGED_TR_DELTA) {
201: snes->reason = SNES_DIVERGED_INNER;
202: PetscFunctionReturn(PETSC_SUCCESS);
203: }
204: PetscCall(SNESGetNPCFunction(snes, F, &fnorm));
205: } else if (snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
206: PetscCall(SNESApplyNPC(snes, X, F, F));
207: PetscCall(SNESGetConvergedReason(snes->npc, &reason));
208: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT && reason != SNES_DIVERGED_TR_DELTA) {
209: snes->reason = SNES_DIVERGED_INNER;
210: PetscFunctionReturn(PETSC_SUCCESS);
211: }
212: }
213: }
215: /* Solve J Y = F, where J is Jacobian matrix */
216: PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre));
217: SNESCheckJacobianDomainerror(snes);
218: PetscCall(KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre));
219: PetscCall(KSPSolve(snes->ksp, F, Y));
220: SNESCheckKSPSolve(snes);
221: PetscCall(KSPGetIterationNumber(snes->ksp, &lits));
222: PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits));
224: if (PetscLogPrintInfo) PetscCall(SNESNEWTONLSCheckResidual_Private(snes, snes->jacobian, F, Y));
226: #if defined(PETSC_USE_INFO)
227: gnorm = fnorm;
228: #endif
229: /* Compute a (scaled) negative update in the line search routine:
230: X <- X - lambda*Y
231: and evaluate F = function(X) (depends on the line search).
232: */
233: PetscCall(SNESLineSearchApply(linesearch, X, F, &fnorm, Y));
234: PetscCall(SNESLineSearchGetReason(linesearch, &lssucceed));
235: PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm));
236: PetscCall(PetscInfo(snes, "fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n", (double)gnorm, (double)fnorm, (double)ynorm, (int)lssucceed));
237: if (snes->reason) break;
238: SNESCheckFunctionNorm(snes, fnorm);
239: if (lssucceed) {
240: if (snes->stol * xnorm > ynorm) {
241: snes->reason = SNES_CONVERGED_SNORM_RELATIVE;
242: PetscFunctionReturn(PETSC_SUCCESS);
243: }
244: if (++snes->numFailures >= snes->maxFailures) {
245: PetscBool ismin;
246: snes->reason = SNES_DIVERGED_LINE_SEARCH;
247: PetscCall(SNESNEWTONLSCheckLocalMin_Private(snes, snes->jacobian, F, fnorm, &ismin));
248: if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
249: if (snes->errorifnotconverged && snes->reason) {
250: PetscViewer monitor;
251: PetscCall(SNESLineSearchGetDefaultMonitor(linesearch, &monitor));
252: PetscCheck(monitor, PetscObjectComm((PetscObject)snes), PETSC_ERR_NOT_CONVERGED, "SNESSolve has not converged due to %s. Suggest running with -snes_linesearch_monitor", SNESConvergedReasons[snes->reason]);
253: SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_NOT_CONVERGED, "SNESSolve has not converged due %s.", SNESConvergedReasons[snes->reason]);
254: }
255: break;
256: }
257: }
258: /* Monitor convergence */
259: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
260: snes->iter = i + 1;
261: snes->norm = fnorm;
262: snes->ynorm = ynorm;
263: snes->xnorm = xnorm;
264: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
265: PetscCall(SNESLogConvergenceHistory(snes, snes->norm, lits));
266: PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
267: /* Test for convergence */
268: PetscUseTypeMethod(snes, converged, snes->iter, xnorm, ynorm, fnorm, &snes->reason, snes->cnvP);
269: if (snes->reason) break;
270: }
271: if (i == maxits) {
272: PetscCall(PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", maxits));
273: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
274: }
275: PetscFunctionReturn(PETSC_SUCCESS);
276: }
278: /*
279: SNESSetUp_NEWTONLS - Sets up the internal data structures for the later use
280: of the SNESNEWTONLS nonlinear solver.
282: Input Parameter:
283: . snes - the SNES context
284: . x - the solution vector
286: Application Interface Routine: SNESSetUp()
288: */
289: PetscErrorCode SNESSetUp_NEWTONLS(SNES snes)
290: {
291: PetscFunctionBegin;
292: PetscCall(SNESSetUpMatrices(snes));
293: if (snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_DEFAULT) snes->functype = SNES_FUNCTION_PRECONDITIONED;
294: PetscFunctionReturn(PETSC_SUCCESS);
295: }
297: PetscErrorCode SNESReset_NEWTONLS(SNES snes)
298: {
299: PetscFunctionBegin;
300: PetscFunctionReturn(PETSC_SUCCESS);
301: }
303: /*
304: SNESDestroy_NEWTONLS - Destroys the private SNES_NEWTONLS context that was created
305: with SNESCreate_NEWTONLS().
307: Input Parameter:
308: . snes - the SNES context
310: Application Interface Routine: SNESDestroy()
311: */
312: PetscErrorCode SNESDestroy_NEWTONLS(SNES snes)
313: {
314: PetscFunctionBegin;
315: PetscCall(SNESReset_NEWTONLS(snes));
316: PetscCall(PetscFree(snes->data));
317: PetscFunctionReturn(PETSC_SUCCESS);
318: }
320: /*
321: SNESView_NEWTONLS - Prints info from the SNESNEWTONLS data structure.
323: Input Parameters:
324: . SNES - the SNES context
325: . viewer - visualization context
327: Application Interface Routine: SNESView()
328: */
329: static PetscErrorCode SNESView_NEWTONLS(SNES snes, PetscViewer viewer)
330: {
331: PetscBool iascii;
333: PetscFunctionBegin;
334: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
335: if (iascii) { }
336: PetscFunctionReturn(PETSC_SUCCESS);
337: }
339: /*
340: SNESSetFromOptions_NEWTONLS - Sets various parameters for the SNESNEWTONLS method.
342: Input Parameter:
343: . snes - the SNES context
345: Application Interface Routine: SNESSetFromOptions()
346: */
347: static PetscErrorCode SNESSetFromOptions_NEWTONLS(SNES snes, PetscOptionItems *PetscOptionsObject)
348: {
349: PetscFunctionBegin;
350: PetscFunctionReturn(PETSC_SUCCESS);
351: }
353: /*MC
354: SNESNEWTONLS - Newton based nonlinear solver that uses a line search
356: Options Database Keys:
357: + -snes_linesearch_type <bt> - bt,basic. Select line search type
358: . -snes_linesearch_order <3> - 2, 3. Selects the order of the line search for bt
359: . -snes_linesearch_norms <true> - Turns on/off computation of the norms for basic linesearch (`SNESLineSearchSetComputeNorms()`)
360: . -snes_linesearch_alpha <alpha> - Sets alpha used in determining if reduction in function norm is sufficient
361: . -snes_linesearch_maxstep <maxstep> - Sets the maximum stepsize the line search will use (if the 2-norm(y) > maxstep then scale y to be y = (maxstep/2-norm(y)) *y)
362: . -snes_linesearch_minlambda <minlambda> - Sets the minimum lambda the line search will tolerate
363: . -snes_linesearch_monitor - print information about progress of line searches
364: - -snes_linesearch_damping - damping factor used for basic line search
366: Note:
367: This is the default nonlinear solver in `SNES`
369: Level: beginner
371: .seealso: `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONTR`, `SNESQN`, `SNESLineSearchSetType()`, `SNESLineSearchSetOrder()`
372: `SNESLineSearchSetPostCheck()`, `SNESLineSearchSetPreCheck()` `SNESLineSearchSetComputeNorms()`, `SNESGetLineSearch()`
373: M*/
374: PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONLS(SNES snes)
375: {
376: SNES_NEWTONLS *neP;
377: SNESLineSearch linesearch;
379: PetscFunctionBegin;
380: snes->ops->setup = SNESSetUp_NEWTONLS;
381: snes->ops->solve = SNESSolve_NEWTONLS;
382: snes->ops->destroy = SNESDestroy_NEWTONLS;
383: snes->ops->setfromoptions = SNESSetFromOptions_NEWTONLS;
384: snes->ops->view = SNESView_NEWTONLS;
385: snes->ops->reset = SNESReset_NEWTONLS;
387: snes->npcside = PC_RIGHT;
388: snes->usesksp = PETSC_TRUE;
389: snes->usesnpc = PETSC_TRUE;
391: PetscCall(SNESGetLineSearch(snes, &linesearch));
392: if (!((PetscObject)linesearch)->type_name) PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT));
394: snes->alwayscomputesfinalresidual = PETSC_TRUE;
396: PetscCall(PetscNew(&neP));
397: snes->data = (void *)neP;
398: PetscFunctionReturn(PETSC_SUCCESS);
399: }