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