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