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