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) {
66: PetscInfo(snes,"||J^T(F-Ax)||/||F-AX|| %14.12e near zero implies inconsistent rhs\n",(double)(a2/a1));
67: }
68: VecDestroy(&W1);
69: VecDestroy(&W2);
70: }
71: return 0;
72: }
74: /* --------------------------------------------------------------------
76: This file implements a truncated Newton method with a line search,
77: for solving a system of nonlinear equations, using the KSP, Vec,
78: and Mat interfaces for linear solvers, vectors, and matrices,
79: respectively.
81: The following basic routines are required for each nonlinear solver:
82: SNESCreate_XXX() - Creates a nonlinear solver context
83: SNESSetFromOptions_XXX() - Sets runtime options
84: SNESSolve_XXX() - Solves the nonlinear system
85: SNESDestroy_XXX() - Destroys the nonlinear solver context
86: The suffix "_XXX" denotes a particular implementation, in this case
87: we use _NEWTONLS (e.g., SNESCreate_NEWTONLS, SNESSolve_NEWTONLS) for solving
88: systems of nonlinear equations with a line search (LS) method.
89: These routines are actually called via the common user interface
90: routines SNESCreate(), SNESSetFromOptions(), SNESSolve(), and
91: SNESDestroy(), so the application code interface remains identical
92: for all nonlinear solvers.
94: Another key routine is:
95: SNESSetUp_XXX() - Prepares for the use of a nonlinear solver
96: by setting data structures and options. The interface routine SNESSetUp()
97: is not usually called directly by the user, but instead is called by
98: SNESSolve() if necessary.
100: Additional basic routines are:
101: SNESView_XXX() - Prints details of runtime options that
102: have actually been used.
103: These are called by application codes via the interface routines
104: SNESView().
106: The various types of solvers (preconditioners, Krylov subspace methods,
107: nonlinear solvers, timesteppers) are all organized similarly, so the
108: above description applies to these categories also.
110: -------------------------------------------------------------------- */
111: /*
112: SNESSolve_NEWTONLS - Solves a nonlinear system with a truncated Newton
113: method with a line search.
115: Input Parameters:
116: . snes - the SNES context
118: Output Parameter:
119: . outits - number of iterations until termination
121: Application Interface Routine: SNESSolve()
123: Notes:
124: This implements essentially a truncated Newton method with a
125: line search. By default a cubic backtracking line search
126: is employed, as described in the text "Numerical Methods for
127: Unconstrained Optimization and Nonlinear Equations" by Dennis
128: and Schnabel.
129: */
130: PetscErrorCode SNESSolve_NEWTONLS(SNES snes)
131: {
132: PetscInt maxits,i,lits;
133: SNESLineSearchReason lssucceed;
134: PetscReal fnorm,gnorm,xnorm,ynorm;
135: Vec Y,X,F;
136: SNESLineSearch linesearch;
137: SNESConvergedReason reason;
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: PetscObjectSAWsTakeAccess((PetscObject)snes);
151: snes->iter = 0;
152: snes->norm = 0.0;
153: PetscObjectSAWsGrantAccess((PetscObject)snes);
154: 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: SNESApplyNPC(snes,X,NULL,F);
159: SNESGetConvergedReason(snes->npc,&reason);
160: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
161: snes->reason = SNES_DIVERGED_INNER;
162: return 0;
163: }
165: VecNormBegin(F,NORM_2,&fnorm);
166: VecNormEnd(F,NORM_2,&fnorm);
167: } else {
168: if (!snes->vec_func_init_set) {
169: SNESComputeFunction(snes,X,F);
170: } else snes->vec_func_init_set = PETSC_FALSE;
171: }
173: VecNorm(F,NORM_2,&fnorm); /* fnorm <- ||F|| */
174: SNESCheckFunctionNorm(snes,fnorm);
175: PetscObjectSAWsTakeAccess((PetscObject)snes);
176: snes->norm = fnorm;
177: PetscObjectSAWsGrantAccess((PetscObject)snes);
178: SNESLogConvergenceHistory(snes,fnorm,0);
179: SNESMonitor(snes,0,fnorm);
181: /* test convergence */
182: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
183: if (snes->reason) return 0;
185: for (i=0; i<maxits; i++) {
187: /* Call general purpose update function */
188: if (snes->ops->update) {
189: (*snes->ops->update)(snes, snes->iter);
190: }
192: /* apply the nonlinear preconditioner */
193: if (snes->npc) {
194: if (snes->npcside== PC_RIGHT) {
195: SNESSetInitialFunction(snes->npc, F);
196: PetscLogEventBegin(SNES_NPCSolve,snes->npc,X,snes->vec_rhs,0);
197: SNESSolve(snes->npc, snes->vec_rhs, X);
198: PetscLogEventEnd(SNES_NPCSolve,snes->npc,X,snes->vec_rhs,0);
199: SNESGetConvergedReason(snes->npc,&reason);
200: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
201: snes->reason = SNES_DIVERGED_INNER;
202: return 0;
203: }
204: SNESGetNPCFunction(snes,F,&fnorm);
205: } else if (snes->npcside== PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
206: SNESApplyNPC(snes,X,F,F);
207: SNESGetConvergedReason(snes->npc,&reason);
208: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
209: snes->reason = SNES_DIVERGED_INNER;
210: return 0;
211: }
212: }
213: }
215: /* Solve J Y = F, where J is Jacobian matrix */
216: SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
217: SNESCheckJacobianDomainerror(snes);
218: KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);
219: KSPSolve(snes->ksp,F,Y);
220: SNESCheckKSPSolve(snes);
221: KSPGetIterationNumber(snes->ksp,&lits);
222: PetscInfo(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);
224: if (PetscLogPrintInfo) {
225: SNESNEWTONLSCheckResidual_Private(snes,snes->jacobian,F,Y);
226: }
228: /* Compute a (scaled) negative update in the line search routine:
229: X <- X - lambda*Y
230: and evaluate F = function(X) (depends on the line search).
231: */
232: gnorm = fnorm;
233: SNESLineSearchApply(linesearch, X, F, &fnorm, Y);
234: SNESLineSearchGetReason(linesearch, &lssucceed);
235: SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm);
236: 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: return 0;
243: }
244: if (++snes->numFailures >= snes->maxFailures) {
245: PetscBool ismin;
246: snes->reason = SNES_DIVERGED_LINE_SEARCH;
247: 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: SNESLineSearchGetDefaultMonitor(linesearch,&monitor);
253: else 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: PetscObjectSAWsTakeAccess((PetscObject)snes);
260: snes->iter = i+1;
261: snes->norm = fnorm;
262: snes->ynorm = ynorm;
263: snes->xnorm = xnorm;
264: PetscObjectSAWsGrantAccess((PetscObject)snes);
265: SNESLogConvergenceHistory(snes,snes->norm,lits);
266: SNESMonitor(snes,snes->iter,snes->norm);
267: /* Test for convergence */
268: (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);
269: if (snes->reason) break;
270: }
271: if (i == maxits) {
272: PetscInfo(snes,"Maximum number of iterations has been reached: %D\n",maxits);
273: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
274: }
275: return 0;
276: }
277: /* -------------------------------------------------------------------------- */
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: 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: */
293: PetscErrorCode SNESSetUp_NEWTONLS(SNES snes)
294: {
295: SNESSetUpMatrices(snes);
296: if (snes->npcside== PC_LEFT && snes->functype == SNES_FUNCTION_DEFAULT) snes->functype = SNES_FUNCTION_PRECONDITIONED;
297: return 0;
298: }
299: /* -------------------------------------------------------------------------- */
301: PetscErrorCode SNESReset_NEWTONLS(SNES snes)
302: {
303: return 0;
304: }
306: /*
307: SNESDestroy_NEWTONLS - Destroys the private SNES_NEWTONLS context that was created
308: with SNESCreate_NEWTONLS().
310: Input Parameter:
311: . snes - the SNES context
313: Application Interface Routine: SNESDestroy()
314: */
315: PetscErrorCode SNESDestroy_NEWTONLS(SNES snes)
316: {
317: SNESReset_NEWTONLS(snes);
318: PetscFree(snes->data);
319: return 0;
320: }
321: /* -------------------------------------------------------------------------- */
323: /*
324: SNESView_NEWTONLS - Prints info from the SNESNEWTONLS data structure.
326: Input Parameters:
327: . SNES - the SNES context
328: . viewer - visualization context
330: Application Interface Routine: SNESView()
331: */
332: static PetscErrorCode SNESView_NEWTONLS(SNES snes,PetscViewer viewer)
333: {
334: PetscBool iascii;
336: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
337: if (iascii) {
338: }
339: return 0;
340: }
342: /* -------------------------------------------------------------------------- */
343: /*
344: SNESSetFromOptions_NEWTONLS - Sets various parameters for the SNESNEWTONLS method.
346: Input Parameter:
347: . snes - the SNES context
349: Application Interface Routine: SNESSetFromOptions()
350: */
351: static PetscErrorCode SNESSetFromOptions_NEWTONLS(PetscOptionItems *PetscOptionsObject,SNES snes)
352: {
353: return 0;
354: }
356: /* -------------------------------------------------------------------------- */
357: /*MC
358: SNESNEWTONLS - Newton based nonlinear solver that uses a line search
360: Options Database:
361: + -snes_linesearch_type <bt> - bt,basic. Select line search type
362: . -snes_linesearch_order <3> - 2, 3. Selects the order of the line search for bt
363: . -snes_linesearch_norms <true> - Turns on/off computation of the norms for basic linesearch (SNESLineSearchSetComputeNorms())
364: . -snes_linesearch_alpha <alpha> - Sets alpha used in determining if reduction in function norm is sufficient
365: . -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)
366: . -snes_linesearch_minlambda <minlambda> - Sets the minimum lambda the line search will tolerate
367: . -snes_linesearch_monitor - print information about progress of line searches
368: - -snes_linesearch_damping - damping factor used for basic line search
370: Notes:
371: This is the default nonlinear solver in SNES
373: Level: beginner
375: .seealso: SNESCreate(), SNES, SNESSetType(), SNESNEWTONTR, SNESQN, SNESLineSearchSetType(), SNESLineSearchSetOrder()
376: SNESLineSearchSetPostCheck(), SNESLineSearchSetPreCheck() SNESLineSearchSetComputeNorms()
378: M*/
379: PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONLS(SNES snes)
380: {
381: SNES_NEWTONLS *neP;
382: SNESLineSearch linesearch;
384: snes->ops->setup = SNESSetUp_NEWTONLS;
385: snes->ops->solve = SNESSolve_NEWTONLS;
386: snes->ops->destroy = SNESDestroy_NEWTONLS;
387: snes->ops->setfromoptions = SNESSetFromOptions_NEWTONLS;
388: snes->ops->view = SNESView_NEWTONLS;
389: snes->ops->reset = SNESReset_NEWTONLS;
391: snes->npcside = PC_RIGHT;
392: snes->usesksp = PETSC_TRUE;
393: snes->usesnpc = PETSC_TRUE;
395: SNESGetLineSearch(snes, &linesearch);
396: if (!((PetscObject)linesearch)->type_name) {
397: SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);
398: }
400: snes->alwayscomputesfinalresidual = PETSC_TRUE;
402: PetscNewLog(snes,&neP);
403: snes->data = (void*)neP;
404: return 0;
405: }