Actual source code: ls.c

petsc-master 2018-06-18
Report Typos and Errors
```
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.

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:     KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);
224:     KSPSolve(snes->ksp,F,Y);
225:     SNESCheckKSPSolve(snes);
226:     KSPGetIterationNumber(snes->ksp,&lits);
227:     snes->linear_its += 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 == SNES_DIVERGED_FUNCTION_COUNT) 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:         break;
256:       }
257:     }
258:     /* Monitor convergence */
259:     PetscObjectSAWsTakeAccess((PetscObject)snes);
260:     snes->iter = i+1;
261:     snes->norm = fnorm;
262:     PetscObjectSAWsGrantAccess((PetscObject)snes);
263:     SNESLogConvergenceHistory(snes,snes->norm,lits);
264:     SNESMonitor(snes,snes->iter,snes->norm);
265:     /* Test for convergence */
266:     (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);
267:     if (snes->reason) break;
268:   }
269:   if (i == maxits) {
270:     PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);
271:     if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
272:   }
273:   return(0);
274: }
275: /* -------------------------------------------------------------------------- */
276: /*
277:    SNESSetUp_NEWTONLS - Sets up the internal data structures for the later use
278:    of the SNESNEWTONLS nonlinear solver.

280:    Input Parameter:
281: .  snes - the SNES context
282: .  x - the solution vector

284:    Application Interface Routine: SNESSetUp()

286:    Notes:
287:    For basic use of the SNES solvers, the user need not explicitly call
288:    SNESSetUp(), since these actions will automatically occur during
289:    the call to SNESSolve().
290:  */
291: PetscErrorCode SNESSetUp_NEWTONLS(SNES snes)
292: {

296:   SNESSetUpMatrices(snes);
297:   if (snes->npcside== PC_LEFT && snes->functype == SNES_FUNCTION_DEFAULT) snes->functype = SNES_FUNCTION_PRECONDITIONED;
298:   return(0);
299: }
300: /* -------------------------------------------------------------------------- */

302: PetscErrorCode SNESReset_NEWTONLS(SNES snes)
303: {
305:   return(0);
306: }

308: /*
309:    SNESDestroy_NEWTONLS - Destroys the private SNES_NEWTONLS context that was created
310:    with SNESCreate_NEWTONLS().

312:    Input Parameter:
313: .  snes - the SNES context

315:    Application Interface Routine: SNESDestroy()
316:  */
317: PetscErrorCode SNESDestroy_NEWTONLS(SNES snes)
318: {

322:   SNESReset_NEWTONLS(snes);
323:   PetscFree(snes->data);
324:   return(0);
325: }
326: /* -------------------------------------------------------------------------- */

328: /*
329:    SNESView_NEWTONLS - Prints info from the SNESNEWTONLS data structure.

331:    Input Parameters:
332: .  SNES - the SNES context
333: .  viewer - visualization context

335:    Application Interface Routine: SNESView()
336: */
337: static PetscErrorCode SNESView_NEWTONLS(SNES snes,PetscViewer viewer)
338: {
340:   PetscBool      iascii;

343:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
344:   if (iascii) {
345:   }
346:   return(0);
347: }

349: /* -------------------------------------------------------------------------- */
350: /*
351:    SNESSetFromOptions_NEWTONLS - Sets various parameters for the SNESNEWTONLS method.

353:    Input Parameter:
354: .  snes - the SNES context

356:    Application Interface Routine: SNESSetFromOptions()
357: */
358: static PetscErrorCode SNESSetFromOptions_NEWTONLS(PetscOptionItems *PetscOptionsObject,SNES snes)
359: {
361:   SNESLineSearch linesearch;

364:   if (!snes->linesearch) {
365:     SNESGetLineSearch(snes, &linesearch);
366:     SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);
367:   }
368:   return(0);
369: }

371: /* -------------------------------------------------------------------------- */
372: /*MC
373:       SNESNEWTONLS - Newton based nonlinear solver that uses a line search

375:    Options Database:
376: +   -snes_linesearch_type <bt> - bt,basic.  Select line search type
377: .   -snes_linesearch_order <3> - 2, 3. Selects the order of the line search for bt
378: .   -snes_linesearch_norms <true> - Turns on/off computation of the norms for basic linesearch (SNESLineSearchSetComputeNorms())
379: .   -snes_linesearch_alpha <alpha> - Sets alpha used in determining if reduction in function norm is sufficient
380: .   -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)
381: .   -snes_linesearch_minlambda <minlambda>  - Sets the minimum lambda the line search will tolerate
382: .   -snes_linesearch_monitor - print information about progress of line searches
383: -   -snes_linesearch_damping - damping factor used for basic line search

385:     Notes:
386:     This is the default nonlinear solver in SNES

388:    Level: beginner

390: .seealso:  SNESCreate(), SNES, SNESSetType(), SNESNEWTONTR, SNESQN, SNESLineSearchSetType(), SNESLineSearchSetOrder()
391:            SNESLineSearchSetPostCheck(), SNESLineSearchSetPreCheck() SNESLineSearchSetComputeNorms()

393: M*/
394: PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONLS(SNES snes)
395: {
397:   SNES_NEWTONLS  *neP;

400:   snes->ops->setup          = SNESSetUp_NEWTONLS;
401:   snes->ops->solve          = SNESSolve_NEWTONLS;
402:   snes->ops->destroy        = SNESDestroy_NEWTONLS;
403:   snes->ops->setfromoptions = SNESSetFromOptions_NEWTONLS;
404:   snes->ops->view           = SNESView_NEWTONLS;
405:   snes->ops->reset          = SNESReset_NEWTONLS;

407:   snes->npcside = PC_RIGHT;
408:   snes->usesksp = PETSC_TRUE;
409:   snes->usesnpc = PETSC_TRUE;

411:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

413:   PetscNewLog(snes,&neP);
414:   snes->data    = (void*)neP;
415:   return(0);
416: }
```