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;
 15:   PetscErrorCode (*objective)(SNES, Vec, PetscReal *, void *);

 17:   PetscFunctionBegin;
 18:   *ismin = PETSC_FALSE;
 19:   PetscCall(SNESGetObjective(snes, &objective, NULL));
 20:   if (!objective) {
 21:     PetscCall(MatHasOperation(A, MATOP_MULT_TRANSPOSE, &hastranspose));
 22:     PetscCall(VecDuplicate(F, &W));
 23:     if (hastranspose) {
 24:       /* Compute || J^T F|| */
 25:       PetscCall(MatMultTranspose(A, F, W));
 26:       PetscCall(VecNorm(W, NORM_2, &a1));
 27:       PetscCall(PetscInfo(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:       PetscCall(VecSetRandom(W, NULL));
 35:       PetscCall(VecNorm(W, NORM_2, &wnorm));
 36:       PetscCall(VecDuplicate(W, &work));
 37:       PetscCall(MatMult(A, W, work));
 38:       PetscCall(VecDot(F, work, &result));
 39:       PetscCall(VecDestroy(&work));
 40:       a1 = PetscAbsScalar(result) / (fnorm * wnorm);
 41:       PetscCall(PetscInfo(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:     PetscCall(VecDestroy(&W));
 45:   }
 46:   PetscFunctionReturn(PETSC_SUCCESS);
 47: }

 49: /*
 50:      Checks if J^T(F - J*X) = 0
 51: */
 52: static PetscErrorCode SNESNEWTONLSCheckResidual_Private(SNES snes, Mat A, Vec F, Vec X)
 53: {
 54:   PetscReal a1, a2;
 55:   PetscBool hastranspose;
 56:   PetscErrorCode (*objective)(SNES, Vec, PetscReal *, void *);

 58:   PetscFunctionBegin;
 59:   PetscCall(MatHasOperation(A, MATOP_MULT_TRANSPOSE, &hastranspose));
 60:   PetscCall(SNESGetObjective(snes, &objective, NULL));
 61:   if (hastranspose && !objective) {
 62:     Vec W1, W2;

 64:     PetscCall(VecDuplicate(F, &W1));
 65:     PetscCall(VecDuplicate(F, &W2));
 66:     PetscCall(MatMult(A, X, W1));
 67:     PetscCall(VecAXPY(W1, -1.0, F));

 69:     /* Compute || J^T W|| */
 70:     PetscCall(MatMultTranspose(A, W1, W2));
 71:     PetscCall(VecNorm(W1, NORM_2, &a1));
 72:     PetscCall(VecNorm(W2, NORM_2, &a2));
 73:     if (a1 != 0.0) PetscCall(PetscInfo(snes, "||J^T(F-Ax)||/||F-AX|| %14.12e near zero implies inconsistent rhs\n", (double)(a2 / a1)));
 74:     PetscCall(VecDestroy(&W1));
 75:     PetscCall(VecDestroy(&W2));
 76:   }
 77:   PetscFunctionReturn(PETSC_SUCCESS);
 78: }

 80: /*
 81:      This file implements a truncated Newton method with a line search,
 82:      for solving a system of nonlinear equations, using the KSP, Vec,
 83:      and Mat interfaces for linear solvers, vectors, and matrices,
 84:      respectively.

 86:      The following basic routines are required for each nonlinear solver:
 87:           SNESCreate_XXX()          - Creates a nonlinear solver context
 88:           SNESSetFromOptions_XXX()  - Sets runtime options
 89:           SNESSolve_XXX()           - Solves the nonlinear system
 90:           SNESDestroy_XXX()         - Destroys the nonlinear solver context
 91:      The suffix "_XXX" denotes a particular implementation, in this case
 92:      we use _NEWTONLS (e.g., SNESCreate_NEWTONLS, SNESSolve_NEWTONLS) for solving
 93:      systems of nonlinear equations with a line search (LS) method.
 94:      These routines are actually called via the common user interface
 95:      routines SNESCreate(), SNESSetFromOptions(), SNESSolve(), and
 96:      SNESDestroy(), so the application code interface remains identical
 97:      for all nonlinear solvers.

 99:      Another key routine is:
100:           SNESSetUp_XXX()           - Prepares for the use of a nonlinear solver
101:      by setting data structures and options.   The interface routine SNESSetUp()
102:      is not usually called directly by the user, but instead is called by
103:      SNESSolve() if necessary.

105:      Additional basic routines are:
106:           SNESView_XXX()            - Prints details of runtime options that
107:                                       have actually been used.
108:      These are called by application codes via the interface routines
109:      SNESView().

111:      The various types of solvers (preconditioners, Krylov subspace methods,
112:      nonlinear solvers, timesteppers) are all organized similarly, so the
113:      above description applies to these categories also.

115: */
116: /*
117:    SNESSolve_NEWTONLS - Solves a nonlinear system with a truncated Newton
118:    method with a line search.

120:    Input Parameter:
121: .  snes - the SNES context

123:    Output Parameter:
124: .  outits - number of iterations until termination

126:    Application Interface Routine: SNESSolve()

128: */
129: PetscErrorCode SNESSolve_NEWTONLS(SNES snes)
130: {
131:   PetscInt             maxits, i, lits;
132:   SNESLineSearchReason lssucceed;
133:   PetscReal            fnorm, xnorm, ynorm;
134:   Vec                  Y, X, F;
135:   SNESLineSearch       linesearch;
136:   SNESConvergedReason  reason;
137: #if defined(PETSC_USE_INFO)
138:   PetscReal gnorm;
139: #endif

141:   PetscFunctionBegin;
142:   PetscCheck(!snes->xl && !snes->xu && !snes->ops->computevariablebounds, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);

144:   snes->numFailures            = 0;
145:   snes->numLinearSolveFailures = 0;
146:   snes->reason                 = SNES_CONVERGED_ITERATING;

148:   maxits = snes->max_its;        /* maximum number of iterations */
149:   X      = snes->vec_sol;        /* solution vector */
150:   F      = snes->vec_func;       /* residual vector */
151:   Y      = snes->vec_sol_update; /* newton step */

153:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
154:   snes->iter = 0;
155:   snes->norm = 0.0;
156:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
157:   PetscCall(SNESGetLineSearch(snes, &linesearch));

159:   /* compute the preconditioned function first in the case of left preconditioning with preconditioned function */
160:   if (snes->npc && snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_PRECONDITIONED) {
161:     PetscCall(SNESApplyNPC(snes, X, NULL, F));
162:     PetscCall(SNESGetConvergedReason(snes->npc, &reason));
163:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT && reason != SNES_DIVERGED_TR_DELTA) {
164:       snes->reason = SNES_DIVERGED_INNER;
165:       PetscFunctionReturn(PETSC_SUCCESS);
166:     }

168:     PetscCall(VecNormBegin(F, NORM_2, &fnorm));
169:     PetscCall(VecNormEnd(F, NORM_2, &fnorm));
170:   } else {
171:     if (!snes->vec_func_init_set) {
172:       PetscCall(SNESComputeFunction(snes, X, F));
173:     } else snes->vec_func_init_set = PETSC_FALSE;
174:   }

176:   PetscCall(VecNorm(F, NORM_2, &fnorm)); /* fnorm <- ||F||  */
177:   SNESCheckFunctionNorm(snes, fnorm);
178:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
179:   snes->norm = fnorm;
180:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
181:   PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0));
182:   PetscCall(SNESMonitor(snes, 0, fnorm));

184:   /* test convergence */
185:   PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
186:   if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);

188:   for (i = 0; i < maxits; i++) {
189:     /* Call general purpose update function */
190:     PetscTryTypeMethod(snes, update, snes->iter);

192:     /* apply the nonlinear preconditioner */
193:     if (snes->npc) {
194:       if (snes->npcside == PC_RIGHT) {
195:         PetscCall(SNESSetInitialFunction(snes->npc, F));
196:         PetscCall(PetscLogEventBegin(SNES_NPCSolve, snes->npc, X, snes->vec_rhs, 0));
197:         PetscCall(SNESSolve(snes->npc, snes->vec_rhs, X));
198:         PetscCall(PetscLogEventEnd(SNES_NPCSolve, snes->npc, X, snes->vec_rhs, 0));
199:         PetscCall(SNESGetConvergedReason(snes->npc, &reason));
200:         if (reason < 0 && reason != SNES_DIVERGED_MAX_IT && reason != SNES_DIVERGED_TR_DELTA) {
201:           snes->reason = SNES_DIVERGED_INNER;
202:           PetscFunctionReturn(PETSC_SUCCESS);
203:         }
204:         PetscCall(SNESGetNPCFunction(snes, F, &fnorm));
205:       } else if (snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
206:         PetscCall(SNESApplyNPC(snes, X, F, F));
207:         PetscCall(SNESGetConvergedReason(snes->npc, &reason));
208:         if (reason < 0 && reason != SNES_DIVERGED_MAX_IT && reason != SNES_DIVERGED_TR_DELTA) {
209:           snes->reason = SNES_DIVERGED_INNER;
210:           PetscFunctionReturn(PETSC_SUCCESS);
211:         }
212:       }
213:     }

215:     /* Solve J Y = F, where J is Jacobian matrix */
216:     PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre));
217:     SNESCheckJacobianDomainerror(snes);
218:     PetscCall(KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre));
219:     PetscCall(KSPSolve(snes->ksp, F, Y));
220:     SNESCheckKSPSolve(snes);
221:     PetscCall(KSPGetIterationNumber(snes->ksp, &lits));
222:     PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits));

224:     if (PetscLogPrintInfo) PetscCall(SNESNEWTONLSCheckResidual_Private(snes, snes->jacobian, F, Y));

226: #if defined(PETSC_USE_INFO)
227:     gnorm = fnorm;
228: #endif
229:     /* Compute a (scaled) negative update in the line search routine:
230:          X <- X - lambda*Y
231:        and evaluate F = function(X) (depends on the line search).
232:     */
233:     PetscCall(SNESLineSearchApply(linesearch, X, F, &fnorm, Y));
234:     PetscCall(SNESLineSearchGetReason(linesearch, &lssucceed));
235:     PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm));
236:     PetscCall(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:         PetscFunctionReturn(PETSC_SUCCESS);
243:       }
244:       if (++snes->numFailures >= snes->maxFailures) {
245:         PetscBool ismin;
246:         snes->reason = SNES_DIVERGED_LINE_SEARCH;
247:         PetscCall(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:           PetscCall(SNESLineSearchGetDefaultMonitor(linesearch, &monitor));
252:           PetscCheck(monitor, PetscObjectComm((PetscObject)snes), PETSC_ERR_NOT_CONVERGED, "SNESSolve has not converged due to %s. Suggest running with -snes_linesearch_monitor", SNESConvergedReasons[snes->reason]);
253:           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:     PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
260:     snes->iter  = i + 1;
261:     snes->norm  = fnorm;
262:     snes->ynorm = ynorm;
263:     snes->xnorm = xnorm;
264:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
265:     PetscCall(SNESLogConvergenceHistory(snes, snes->norm, lits));
266:     PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
267:     /* Test for convergence */
268:     PetscUseTypeMethod(snes, converged, snes->iter, xnorm, ynorm, fnorm, &snes->reason, snes->cnvP);
269:     if (snes->reason) break;
270:   }
271:   if (i == maxits) {
272:     PetscCall(PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", maxits));
273:     if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
274:   }
275:   PetscFunctionReturn(PETSC_SUCCESS);
276: }

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:  */
289: PetscErrorCode SNESSetUp_NEWTONLS(SNES snes)
290: {
291:   PetscFunctionBegin;
292:   PetscCall(SNESSetUpMatrices(snes));
293:   if (snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_DEFAULT) snes->functype = SNES_FUNCTION_PRECONDITIONED;
294:   PetscFunctionReturn(PETSC_SUCCESS);
295: }

297: PetscErrorCode SNESReset_NEWTONLS(SNES snes)
298: {
299:   PetscFunctionBegin;
300:   PetscFunctionReturn(PETSC_SUCCESS);
301: }

303: /*
304:    SNESDestroy_NEWTONLS - Destroys the private SNES_NEWTONLS context that was created
305:    with SNESCreate_NEWTONLS().

307:    Input Parameter:
308: .  snes - the SNES context

310:    Application Interface Routine: SNESDestroy()
311:  */
312: PetscErrorCode SNESDestroy_NEWTONLS(SNES snes)
313: {
314:   PetscFunctionBegin;
315:   PetscCall(SNESReset_NEWTONLS(snes));
316:   PetscCall(PetscFree(snes->data));
317:   PetscFunctionReturn(PETSC_SUCCESS);
318: }

320: /*
321:    SNESView_NEWTONLS - Prints info from the SNESNEWTONLS data structure.

323:    Input Parameters:
324: .  SNES - the SNES context
325: .  viewer - visualization context

327:    Application Interface Routine: SNESView()
328: */
329: static PetscErrorCode SNESView_NEWTONLS(SNES snes, PetscViewer viewer)
330: {
331:   PetscBool iascii;

333:   PetscFunctionBegin;
334:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
335:   if (iascii) { }
336:   PetscFunctionReturn(PETSC_SUCCESS);
337: }

339: /*
340:    SNESSetFromOptions_NEWTONLS - Sets various parameters for the SNESNEWTONLS method.

342:    Input Parameter:
343: .  snes - the SNES context

345:    Application Interface Routine: SNESSetFromOptions()
346: */
347: static PetscErrorCode SNESSetFromOptions_NEWTONLS(SNES snes, PetscOptionItems *PetscOptionsObject)
348: {
349:   PetscFunctionBegin;
350:   PetscFunctionReturn(PETSC_SUCCESS);
351: }

353: /*MC
354:       SNESNEWTONLS - Newton based nonlinear solver that uses a line search

356:    Options Database Keys:
357: +   -snes_linesearch_type <bt> - bt,basic.  Select line search type
358: .   -snes_linesearch_order <3> - 2, 3. Selects the order of the line search for bt
359: .   -snes_linesearch_norms <true> - Turns on/off computation of the norms for basic linesearch (`SNESLineSearchSetComputeNorms()`)
360: .   -snes_linesearch_alpha <alpha> - Sets alpha used in determining if reduction in function norm is sufficient
361: .   -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)
362: .   -snes_linesearch_minlambda <minlambda>  - Sets the minimum lambda the line search will tolerate
363: .   -snes_linesearch_monitor - print information about progress of line searches
364: -   -snes_linesearch_damping - damping factor used for basic line search

366:     Note:
367:     This is the default nonlinear solver in `SNES`

369:    Level: beginner

371: .seealso: `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONTR`, `SNESQN`, `SNESLineSearchSetType()`, `SNESLineSearchSetOrder()`
372:           `SNESLineSearchSetPostCheck()`, `SNESLineSearchSetPreCheck()` `SNESLineSearchSetComputeNorms()`, `SNESGetLineSearch()`
373: M*/
374: PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONLS(SNES snes)
375: {
376:   SNES_NEWTONLS *neP;
377:   SNESLineSearch linesearch;

379:   PetscFunctionBegin;
380:   snes->ops->setup          = SNESSetUp_NEWTONLS;
381:   snes->ops->solve          = SNESSolve_NEWTONLS;
382:   snes->ops->destroy        = SNESDestroy_NEWTONLS;
383:   snes->ops->setfromoptions = SNESSetFromOptions_NEWTONLS;
384:   snes->ops->view           = SNESView_NEWTONLS;
385:   snes->ops->reset          = SNESReset_NEWTONLS;

387:   snes->npcside = PC_RIGHT;
388:   snes->usesksp = PETSC_TRUE;
389:   snes->usesnpc = PETSC_TRUE;

391:   PetscCall(SNESGetLineSearch(snes, &linesearch));
392:   if (!((PetscObject)linesearch)->type_name) PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT));

394:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

396:   PetscCall(PetscNew(&neP));
397:   snes->data = (void *)neP;
398:   PetscFunctionReturn(PETSC_SUCCESS);
399: }