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;
14: PetscBool hastranspose;
15: Vec W;
18: *ismin = PETSC_FALSE;
19: MatHasOperation(A,MATOP_MULT_TRANSPOSE,&hastranspose);
20: VecDuplicate(F,&W);
21: if (hastranspose) {
22: /* Compute || J^T F|| */
23: MatMultTranspose(A,F,W);
24: VecNorm(W,NORM_2,&a1);
25: PetscInfo1(snes,"|| J^T F|| %14.12e near zero implies found a local minimum\n",(double)(a1/fnorm));
26: if (a1/fnorm < 1.e-4) *ismin = PETSC_TRUE;
27: } else {
28: Vec work;
29: PetscScalar result;
30: PetscReal wnorm;
32: VecSetRandom(W,NULL);
33: VecNorm(W,NORM_2,&wnorm);
34: VecDuplicate(W,&work);
35: MatMult(A,W,work);
36: VecDot(F,work,&result);
37: VecDestroy(&work);
38: a1 = PetscAbsScalar(result)/(fnorm*wnorm);
39: PetscInfo1(snes,"(F^T J random)/(|| F ||*||random|| %14.12e near zero implies found a local minimum\n",(double)a1);
40: if (a1 < 1.e-4) *ismin = PETSC_TRUE;
41: }
42: VecDestroy(&W);
43: return(0);
44: }
46: /*
47: Checks if J^T(F - J*X) = 0
48: */
49: static PetscErrorCode SNESNEWTONLSCheckResidual_Private(SNES snes,Mat A,Vec F,Vec X)
50: {
51: PetscReal a1,a2;
53: PetscBool hastranspose;
56: MatHasOperation(A,MATOP_MULT_TRANSPOSE,&hastranspose);
57: if (hastranspose) {
58: Vec W1,W2;
60: VecDuplicate(F,&W1);
61: VecDuplicate(F,&W2);
62: MatMult(A,X,W1);
63: VecAXPY(W1,-1.0,F);
65: /* Compute || J^T W|| */
66: MatMultTranspose(A,W1,W2);
67: VecNorm(W1,NORM_2,&a1);
68: VecNorm(W2,NORM_2,&a2);
69: if (a1 != 0.0) {
70: PetscInfo1(snes,"||J^T(F-Ax)||/||F-AX|| %14.12e near zero implies inconsistent rhs\n",(double)(a2/a1));
71: }
72: VecDestroy(&W1);
73: VecDestroy(&W2);
74: }
75: return(0);
76: }
78: /* --------------------------------------------------------------------
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.
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.
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.
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().
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.
114: -------------------------------------------------------------------- */
115: /*
116: SNESSolve_NEWTONLS - Solves a nonlinear system with a truncated Newton
117: method with a line search.
119: Input Parameters:
120: . snes - the SNES context
122: Output Parameter:
123: . outits - number of iterations until termination
125: Application Interface Routine: SNESSolve()
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: */
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;
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);
147: snes->numFailures = 0;
148: snes->numLinearSolveFailures = 0;
149: snes->reason = SNES_CONVERGED_ITERATING;
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 */
156: PetscObjectSAWsTakeAccess((PetscObject)snes);
157: snes->iter = 0;
158: snes->norm = 0.0;
159: PetscObjectSAWsGrantAccess((PetscObject)snes);
160: SNESGetLineSearch(snes, &linesearch);
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: SNESApplyNPC(snes,X,NULL,F);
165: SNESGetConvergedReason(snes->npc,&reason);
166: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
167: snes->reason = SNES_DIVERGED_INNER;
168: return(0);
169: }
171: VecNormBegin(F,NORM_2,&fnorm);
172: VecNormEnd(F,NORM_2,&fnorm);
173: } else {
174: if (!snes->vec_func_init_set) {
175: SNESComputeFunction(snes,X,F);
176: } else snes->vec_func_init_set = PETSC_FALSE;
177: }
179: VecNorm(F,NORM_2,&fnorm); /* fnorm <- ||F|| */
180: SNESCheckFunctionNorm(snes,fnorm);
181: PetscObjectSAWsTakeAccess((PetscObject)snes);
182: snes->norm = fnorm;
183: PetscObjectSAWsGrantAccess((PetscObject)snes);
184: SNESLogConvergenceHistory(snes,fnorm,0);
185: SNESMonitor(snes,0,fnorm);
187: /* test convergence */
188: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
189: if (snes->reason) return(0);
191: for (i=0; i<maxits; i++) {
193: /* Call general purpose update function */
194: if (snes->ops->update) {
195: (*snes->ops->update)(snes, snes->iter);
196: }
198: /* apply the nonlinear preconditioner */
199: if (snes->npc) {
200: if (snes->npcside== PC_RIGHT) {
201: SNESSetInitialFunction(snes->npc, F);
202: PetscLogEventBegin(SNES_NPCSolve,snes->npc,X,snes->vec_rhs,0);
203: SNESSolve(snes->npc, snes->vec_rhs, X);
204: PetscLogEventEnd(SNES_NPCSolve,snes->npc,X,snes->vec_rhs,0);
205: SNESGetConvergedReason(snes->npc,&reason);
206: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
207: snes->reason = SNES_DIVERGED_INNER;
208: return(0);
209: }
210: SNESGetNPCFunction(snes,F,&fnorm);
211: } else if (snes->npcside== PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
212: SNESApplyNPC(snes,X,F,F);
213: SNESGetConvergedReason(snes->npc,&reason);
214: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
215: snes->reason = SNES_DIVERGED_INNER;
216: return(0);
217: }
218: }
219: }
221: /* Solve J Y = F, where J is Jacobian matrix */
222: SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
223: SNESCheckJacobianDomainerror(snes);
224: KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);
225: KSPSolve(snes->ksp,F,Y);
226: SNESCheckKSPSolve(snes);
227: KSPGetIterationNumber(snes->ksp,&lits);
228: PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);
230: if (PetscLogPrintInfo) {
231: SNESNEWTONLSCheckResidual_Private(snes,snes->jacobian,F,Y);
232: }
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: SNESLineSearchApply(linesearch, X, F, &fnorm, Y);
240: SNESLineSearchGetReason(linesearch, &lssucceed);
241: SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm);
242: PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",(double)gnorm,(double)fnorm,(double)ynorm,(int)lssucceed);
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: return(0);
249: }
250: if (++snes->numFailures >= snes->maxFailures) {
251: PetscBool ismin;
252: snes->reason = SNES_DIVERGED_LINE_SEARCH;
253: SNESNEWTONLSCheckLocalMin_Private(snes,snes->jacobian,F,fnorm,&ismin);
254: if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
255: if (snes->errorifnotconverged && snes->reason) {
256: PetscViewer monitor;
257: SNESLineSearchGetDefaultMonitor(linesearch,&monitor);
258: if (!monitor) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged due to %s. Suggest running with -snes_linesearch_monitor",SNESConvergedReasons[snes->reason]);
259: else SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged due %s.",SNESConvergedReasons[snes->reason]);
260: }
261: break;
262: }
263: }
264: /* Monitor convergence */
265: PetscObjectSAWsTakeAccess((PetscObject)snes);
266: snes->iter = i+1;
267: snes->norm = fnorm;
268: snes->ynorm = ynorm;
269: snes->xnorm = xnorm;
270: PetscObjectSAWsGrantAccess((PetscObject)snes);
271: SNESLogConvergenceHistory(snes,snes->norm,lits);
272: SNESMonitor(snes,snes->iter,snes->norm);
273: /* Test for convergence */
274: (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);
275: if (snes->reason) break;
276: }
277: if (i == maxits) {
278: PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);
279: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
280: }
281: return(0);
282: }
283: /* -------------------------------------------------------------------------- */
284: /*
285: SNESSetUp_NEWTONLS - Sets up the internal data structures for the later use
286: of the SNESNEWTONLS nonlinear solver.
288: Input Parameter:
289: . snes - the SNES context
290: . x - the solution vector
292: Application Interface Routine: SNESSetUp()
294: Notes:
295: For basic use of the SNES solvers, the user need not explicitly call
296: SNESSetUp(), since these actions will automatically occur during
297: the call to SNESSolve().
298: */
299: PetscErrorCode SNESSetUp_NEWTONLS(SNES snes)
300: {
304: SNESSetUpMatrices(snes);
305: if (snes->npcside== PC_LEFT && snes->functype == SNES_FUNCTION_DEFAULT) snes->functype = SNES_FUNCTION_PRECONDITIONED;
306: return(0);
307: }
308: /* -------------------------------------------------------------------------- */
310: PetscErrorCode SNESReset_NEWTONLS(SNES snes)
311: {
313: return(0);
314: }
316: /*
317: SNESDestroy_NEWTONLS - Destroys the private SNES_NEWTONLS context that was created
318: with SNESCreate_NEWTONLS().
320: Input Parameter:
321: . snes - the SNES context
323: Application Interface Routine: SNESDestroy()
324: */
325: PetscErrorCode SNESDestroy_NEWTONLS(SNES snes)
326: {
330: SNESReset_NEWTONLS(snes);
331: PetscFree(snes->data);
332: return(0);
333: }
334: /* -------------------------------------------------------------------------- */
336: /*
337: SNESView_NEWTONLS - Prints info from the SNESNEWTONLS data structure.
339: Input Parameters:
340: . SNES - the SNES context
341: . viewer - visualization context
343: Application Interface Routine: SNESView()
344: */
345: static PetscErrorCode SNESView_NEWTONLS(SNES snes,PetscViewer viewer)
346: {
348: PetscBool iascii;
351: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
352: if (iascii) {
353: }
354: return(0);
355: }
357: /* -------------------------------------------------------------------------- */
358: /*
359: SNESSetFromOptions_NEWTONLS - Sets various parameters for the SNESNEWTONLS method.
361: Input Parameter:
362: . snes - the SNES context
364: Application Interface Routine: SNESSetFromOptions()
365: */
366: static PetscErrorCode SNESSetFromOptions_NEWTONLS(PetscOptionItems *PetscOptionsObject,SNES snes)
367: {
369: return(0);
370: }
372: /* -------------------------------------------------------------------------- */
373: /*MC
374: SNESNEWTONLS - Newton based nonlinear solver that uses a line search
376: Options Database:
377: + -snes_linesearch_type <bt> - bt,basic. Select line search type
378: . -snes_linesearch_order <3> - 2, 3. Selects the order of the line search for bt
379: . -snes_linesearch_norms <true> - Turns on/off computation of the norms for basic linesearch (SNESLineSearchSetComputeNorms())
380: . -snes_linesearch_alpha <alpha> - Sets alpha used in determining if reduction in function norm is sufficient
381: . -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)
382: . -snes_linesearch_minlambda <minlambda> - Sets the minimum lambda the line search will tolerate
383: . -snes_linesearch_monitor - print information about progress of line searches
384: - -snes_linesearch_damping - damping factor used for basic line search
386: Notes:
387: This is the default nonlinear solver in SNES
389: Level: beginner
391: .seealso: SNESCreate(), SNES, SNESSetType(), SNESNEWTONTR, SNESQN, SNESLineSearchSetType(), SNESLineSearchSetOrder()
392: SNESLineSearchSetPostCheck(), SNESLineSearchSetPreCheck() SNESLineSearchSetComputeNorms()
394: M*/
395: PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONLS(SNES snes)
396: {
398: SNES_NEWTONLS *neP;
399: SNESLineSearch linesearch;
402: snes->ops->setup = SNESSetUp_NEWTONLS;
403: snes->ops->solve = SNESSolve_NEWTONLS;
404: snes->ops->destroy = SNESDestroy_NEWTONLS;
405: snes->ops->setfromoptions = SNESSetFromOptions_NEWTONLS;
406: snes->ops->view = SNESView_NEWTONLS;
407: snes->ops->reset = SNESReset_NEWTONLS;
409: snes->npcside = PC_RIGHT;
410: snes->usesksp = PETSC_TRUE;
411: snes->usesnpc = PETSC_TRUE;
413: SNESGetLineSearch(snes, &linesearch);
414: if (!((PetscObject)linesearch)->type_name) {
415: SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);
416: }
418: snes->alwayscomputesfinalresidual = PETSC_TRUE;
420: PetscNewLog(snes,&neP);
421: snes->data = (void*)neP;
422: return(0);
423: }