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