1 #include <../src/snes/impls/ls/lsimpl.h>
2
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.
8
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.
21
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.
27
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().
33
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.
37
38 */
39
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 */
SNESNEWTONLSCheckLocalMin_Private(SNES snes,Mat A,Vec F,PetscReal fnorm,PetscBool * ismin)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;
52
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;
69
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 }
84
85 /*
86 Checks if J^T(F - J*X) = 0
87 */
SNESNEWTONLSCheckResidual_Private(SNES snes,Mat A,Vec F,Vec X)88 static PetscErrorCode SNESNEWTONLSCheckResidual_Private(SNES snes, Mat A, Vec F, Vec X)
89 {
90 PetscReal a1, a2;
91 PetscBool hastranspose;
92 SNESObjectiveFn *objective;
93
94 PetscFunctionBegin;
95 PetscCall(MatHasOperation(A, MATOP_MULT_TRANSPOSE, &hastranspose));
96 PetscCall(SNESGetObjective(snes, &objective, NULL));
97 if (hastranspose && !objective) {
98 Vec W1, W2;
99
100 PetscCall(VecDuplicate(F, &W1));
101 PetscCall(VecDuplicate(F, &W2));
102 PetscCall(MatMult(A, X, W1));
103 PetscCall(VecAXPY(W1, -1.0, F));
104
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 }
115
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.
120
121 Input Parameter:
122 . snes - the SNES context
123
124 */
SNESSolve_NEWTONLS(SNES snes)125 static PetscErrorCode SNESSolve_NEWTONLS(SNES snes)
126 {
127 PetscInt maxits, i, lits;
128 SNESLineSearchReason lsreason;
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
137
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);
140
141 snes->numFailures = 0;
142 snes->numLinearSolveFailures = 0;
143 snes->reason = SNES_CONVERGED_ITERATING;
144
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 */
149
150 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
151 snes->iter = 0;
152 snes->norm = 0.0;
153 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
154 PetscCall(SNESGetLineSearch(snes, &linesearch));
155
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 }
164
165 PetscCall(VecNormBegin(F, NORM_2, &fnorm));
166 PetscCall(VecNormEnd(F, NORM_2, &fnorm));
167 } else {
168 if (!snes->vec_func_init_set) PetscCall(SNESComputeFunction(snes, X, F));
169 else snes->vec_func_init_set = PETSC_FALSE;
170 }
171
172 PetscCall(VecNorm(F, NORM_2, &fnorm)); /* fnorm <- ||F|| */
173 SNESCheckFunctionDomainError(snes, fnorm);
174 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
175 snes->norm = fnorm;
176 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
177 PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0));
178
179 /* test convergence */
180 PetscCall(SNESConverged(snes, 0, 0.0, 0.0, fnorm));
181 PetscCall(SNESMonitor(snes, 0, fnorm));
182 if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
183
184 /* hook state vector to BFGS preconditioner */
185 PetscCall(KSPGetPC(snes->ksp, &pc));
186 PetscCall(PCLMVMSetUpdateVec(pc, X));
187
188 for (i = 0; i < maxits; i++) {
189 /* Call general purpose update function */
190 PetscTryTypeMethod(snes, update, snes->iter);
191 PetscCall(VecNorm(snes->vec_func, NORM_2, &fnorm)); /* no-op unless update() function changed f() */
192
193 /* apply the nonlinear preconditioner */
194 if (snes->npc) {
195 if (snes->npcside == PC_RIGHT) {
196 PetscCall(SNESSetInitialFunction(snes->npc, F));
197 PetscCall(PetscLogEventBegin(SNES_NPCSolve, snes->npc, X, snes->vec_rhs, 0));
198 PetscCall(SNESSolve(snes->npc, snes->vec_rhs, X));
199 PetscCall(PetscLogEventEnd(SNES_NPCSolve, snes->npc, X, snes->vec_rhs, 0));
200 PetscCall(SNESGetConvergedReason(snes->npc, &reason));
201 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT && reason != SNES_DIVERGED_TR_DELTA) {
202 PetscCall(SNESSetConvergedReason(snes, SNES_DIVERGED_INNER));
203 PetscFunctionReturn(PETSC_SUCCESS);
204 }
205 PetscCall(SNESGetNPCFunction(snes, F, &fnorm));
206 } else if (snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
207 PetscCall(SNESApplyNPC(snes, X, F, F));
208 PetscCall(SNESGetConvergedReason(snes->npc, &reason));
209 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT && reason != SNES_DIVERGED_TR_DELTA) {
210 PetscCall(SNESSetConvergedReason(snes, SNES_DIVERGED_INNER));
211 PetscFunctionReturn(PETSC_SUCCESS);
212 }
213 }
214 }
215
216 /* Solve J Y = F, where J is Jacobian matrix */
217 PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre));
218 SNESCheckJacobianDomainError(snes);
219 PetscCall(KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre));
220 PetscCall(KSPSolve(snes->ksp, F, Y));
221 SNESCheckKSPSolve(snes);
222 PetscCall(KSPGetIterationNumber(snes->ksp, &lits));
223 PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits));
224
225 if (PetscLogPrintInfo) PetscCall(SNESNEWTONLSCheckResidual_Private(snes, snes->jacobian, F, Y));
226
227 #if defined(PETSC_USE_INFO)
228 gnorm = fnorm;
229 #endif
230 /* Compute a (scaled) negative update in the line search routine:
231 X <- X - lambda*Y
232 and evaluate F = function(X) (depends on the line search).
233 */
234 PetscCall(SNESLineSearchApply(linesearch, X, F, &fnorm, Y));
235 if (snes->reason) break;
236 PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm));
237 PetscCall(PetscInfo(snes, "fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lsreason=%d\n", (double)gnorm, (double)fnorm, (double)ynorm, (int)lsreason));
238 SNESCheckFunctionDomainError(snes, fnorm);
239 PetscCall(SNESLineSearchGetReason(linesearch, &lsreason));
240 if (lsreason) {
241 if (snes->stol * xnorm > ynorm) {
242 snes->reason = SNES_CONVERGED_SNORM_RELATIVE;
243 break;
244 } else if (lsreason == SNES_LINESEARCH_FAILED_FUNCTION_DOMAIN) {
245 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
246 break;
247 } else if (lsreason == SNES_LINESEARCH_FAILED_NANORINF) {
248 snes->reason = SNES_DIVERGED_FUNCTION_NANORINF;
249 break;
250 } else if (lsreason == SNES_LINESEARCH_FAILED_OBJECTIVE_DOMAIN) {
251 snes->reason = SNES_DIVERGED_OBJECTIVE_DOMAIN;
252 break;
253 } else if (lsreason == SNES_LINESEARCH_FAILED_JACOBIAN_DOMAIN) {
254 snes->reason = SNES_DIVERGED_JACOBIAN_DOMAIN;
255 break;
256 } else if (++snes->numFailures >= snes->maxFailures) {
257 PetscBool ismin;
258
259 snes->reason = SNES_DIVERGED_LINE_SEARCH;
260 PetscCall(SNESNEWTONLSCheckLocalMin_Private(snes, snes->jacobian, F, fnorm, &ismin));
261 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
262 break;
263 }
264 }
265 /* Monitor convergence */
266 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
267 snes->iter = i + 1;
268 snes->norm = fnorm;
269 snes->ynorm = ynorm;
270 snes->xnorm = xnorm;
271 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
272 PetscCall(SNESLogConvergenceHistory(snes, snes->norm, lits));
273 /* Test for convergence */
274 PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, fnorm));
275 PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
276 if (snes->reason) break;
277 }
278 PetscFunctionReturn(PETSC_SUCCESS);
279 }
280
281 /*
282 SNESSetUp_NEWTONLS - Sets up the internal data structures for the later use
283 of the SNESNEWTONLS nonlinear solver.
284
285 Input Parameter:
286 . snes - the SNES context
287 . x - the solution vector
288
289 Application Interface Routine: SNESSetUp()
290
291 */
SNESSetUp_NEWTONLS(SNES snes)292 static PetscErrorCode SNESSetUp_NEWTONLS(SNES snes)
293 {
294 PetscFunctionBegin;
295 PetscCall(SNESSetUpMatrices(snes));
296 if (snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_DEFAULT) snes->functype = SNES_FUNCTION_PRECONDITIONED;
297 PetscFunctionReturn(PETSC_SUCCESS);
298 }
299
300 /*
301 SNESDestroy_NEWTONLS - Destroys the private SNES_NEWTONLS context that was created
302 with SNESCreate_NEWTONLS().
303
304 Input Parameter:
305 . snes - the SNES context
306
307 Application Interface Routine: SNESDestroy()
308 */
SNESDestroy_NEWTONLS(SNES snes)309 static PetscErrorCode SNESDestroy_NEWTONLS(SNES snes)
310 {
311 PetscFunctionBegin;
312 PetscCall(PetscFree(snes->data));
313 PetscFunctionReturn(PETSC_SUCCESS);
314 }
315
316 /*MC
317 SNESNEWTONLS - Newton based nonlinear solver that uses a line search
318
319 Options Database Keys:
320 + -snes_linesearch_type <bt> - basic (or equivalently none), bt, l2, cp, nleqerr, shell. Select line search type, see `SNESLineSearchSetType()`
321 . -snes_linesearch_order <3> - 2, 3. Selects the order of the line search for bt, see `SNESLineSearchSetOrder()`
322 . -snes_linesearch_norms <true> - Turns on/off computation of the norms for basic linesearch (`SNESLineSearchSetComputeNorms()`)
323 . -snes_linesearch_alpha <alpha> - Sets alpha used in determining if reduction in function norm is sufficient
324 . -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)
325 . -snes_linesearch_minlambda <minlambda> - Sets the minimum lambda the line search will tolerate
326 . -snes_linesearch_monitor - print information about the progress of line searches
327 - -snes_linesearch_damping - damping factor used for basic line search
328
329 Level: beginner
330
331 Note:
332 This is the default nonlinear solver in `SNES`
333
334 .seealso: [](ch_snes), `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONTR`, `SNESQN`, `SNESLineSearchSetType()`, `SNESLineSearchSetOrder()`
335 `SNESLineSearchSetPostCheck()`, `SNESLineSearchSetPreCheck()` `SNESLineSearchSetComputeNorms()`, `SNESGetLineSearch()`, `SNESLineSearchSetType()`
336 M*/
SNESCreate_NEWTONLS(SNES snes)337 PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONLS(SNES snes)
338 {
339 SNES_NEWTONLS *neP;
340 SNESLineSearch linesearch;
341
342 PetscFunctionBegin;
343 snes->ops->setup = SNESSetUp_NEWTONLS;
344 snes->ops->solve = SNESSolve_NEWTONLS;
345 snes->ops->destroy = SNESDestroy_NEWTONLS;
346
347 snes->npcside = PC_RIGHT;
348 snes->usesksp = PETSC_TRUE;
349 snes->usesnpc = PETSC_TRUE;
350
351 PetscCall(SNESGetLineSearch(snes, &linesearch));
352 if (!((PetscObject)linesearch)->type_name) PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT));
353
354 snes->alwayscomputesfinalresidual = PETSC_TRUE;
355
356 PetscCall(SNESParametersInitialize(snes));
357
358 PetscCall(PetscNew(&neP));
359 snes->data = (void *)neP;
360 PetscFunctionReturn(PETSC_SUCCESS);
361 }
362