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) PetscInfo(snes, "||J^T(F-Ax)||/||F-AX|| %14.12e near zero implies inconsistent rhs\n", (double)(a2 / a1));
 66:     VecDestroy(&W1);
 67:     VecDestroy(&W2);
 68:   }
 69:   return 0;
 70: }

 72: /*
 73:      This file implements a truncated Newton method with a line search,
 74:      for solving a system of nonlinear equations, using the KSP, Vec,
 75:      and Mat interfaces for linear solvers, vectors, and matrices,
 76:      respectively.

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

 91:      Another key routine is:
 92:           SNESSetUp_XXX()           - Prepares for the use of a nonlinear solver
 93:      by setting data structures and options.   The interface routine SNESSetUp()
 94:      is not usually called directly by the user, but instead is called by
 95:      SNESSolve() if necessary.

 97:      Additional basic routines are:
 98:           SNESView_XXX()            - Prints details of runtime options that
 99:                                       have actually been used.
100:      These are called by application codes via the interface routines
101:      SNESView().

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

107: */
108: /*
109:    SNESSolve_NEWTONLS - Solves a nonlinear system with a truncated Newton
110:    method with a line search.

112:    Input Parameters:
113: .  snes - the SNES context

115:    Output Parameter:
116: .  outits - number of iterations until termination

118:    Application Interface Routine: SNESSolve()

120: */
121: PetscErrorCode SNESSolve_NEWTONLS(SNES snes)
122: {
123:   PetscInt             maxits, i, lits;
124:   SNESLineSearchReason lssucceed;
125:   PetscReal            fnorm, gnorm, xnorm, ynorm;
126:   Vec                  Y, X, F;
127:   SNESLineSearch       linesearch;
128:   SNESConvergedReason  reason;


132:   snes->numFailures            = 0;
133:   snes->numLinearSolveFailures = 0;
134:   snes->reason                 = SNES_CONVERGED_ITERATING;

136:   maxits = snes->max_its;        /* maximum number of iterations */
137:   X      = snes->vec_sol;        /* solution vector */
138:   F      = snes->vec_func;       /* residual vector */
139:   Y      = snes->vec_sol_update; /* newton step */

141:   PetscObjectSAWsTakeAccess((PetscObject)snes);
142:   snes->iter = 0;
143:   snes->norm = 0.0;
144:   PetscObjectSAWsGrantAccess((PetscObject)snes);
145:   SNESGetLineSearch(snes, &linesearch);

147:   /* compute the preconditioned function first in the case of left preconditioning with preconditioned function */
148:   if (snes->npc && snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_PRECONDITIONED) {
149:     SNESApplyNPC(snes, X, NULL, F);
150:     SNESGetConvergedReason(snes->npc, &reason);
151:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
152:       snes->reason = SNES_DIVERGED_INNER;
153:       return 0;
154:     }

156:     VecNormBegin(F, NORM_2, &fnorm);
157:     VecNormEnd(F, NORM_2, &fnorm);
158:   } else {
159:     if (!snes->vec_func_init_set) {
160:       SNESComputeFunction(snes, X, F);
161:     } else snes->vec_func_init_set = PETSC_FALSE;
162:   }

164:   VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
165:   SNESCheckFunctionNorm(snes, fnorm);
166:   PetscObjectSAWsTakeAccess((PetscObject)snes);
167:   snes->norm = fnorm;
168:   PetscObjectSAWsGrantAccess((PetscObject)snes);
169:   SNESLogConvergenceHistory(snes, fnorm, 0);
170:   SNESMonitor(snes, 0, fnorm);

172:   /* test convergence */
173:   PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
174:   if (snes->reason) return 0;

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

180:     /* apply the nonlinear preconditioner */
181:     if (snes->npc) {
182:       if (snes->npcside == PC_RIGHT) {
183:         SNESSetInitialFunction(snes->npc, F);
184:         PetscLogEventBegin(SNES_NPCSolve, snes->npc, X, snes->vec_rhs, 0);
185:         SNESSolve(snes->npc, snes->vec_rhs, X);
186:         PetscLogEventEnd(SNES_NPCSolve, snes->npc, X, snes->vec_rhs, 0);
187:         SNESGetConvergedReason(snes->npc, &reason);
188:         if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
189:           snes->reason = SNES_DIVERGED_INNER;
190:           return 0;
191:         }
192:         SNESGetNPCFunction(snes, F, &fnorm);
193:       } else if (snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
194:         SNESApplyNPC(snes, X, F, F);
195:         SNESGetConvergedReason(snes->npc, &reason);
196:         if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
197:           snes->reason = SNES_DIVERGED_INNER;
198:           return 0;
199:         }
200:       }
201:     }

203:     /* Solve J Y = F, where J is Jacobian matrix */
204:     SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre);
205:     SNESCheckJacobianDomainerror(snes);
206:     KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre);
207:     KSPSolve(snes->ksp, F, Y);
208:     SNESCheckKSPSolve(snes);
209:     KSPGetIterationNumber(snes->ksp, &lits);
210:     PetscInfo(snes, "iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits);

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

214:     /* Compute a (scaled) negative update in the line search routine:
215:          X <- X - lambda*Y
216:        and evaluate F = function(X) (depends on the line search).
217:     */
218:     gnorm = fnorm;
219:     SNESLineSearchApply(linesearch, X, F, &fnorm, Y);
220:     SNESLineSearchGetReason(linesearch, &lssucceed);
221:     SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm);
222:     PetscInfo(snes, "fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n", (double)gnorm, (double)fnorm, (double)ynorm, (int)lssucceed);
223:     if (snes->reason) break;
224:     SNESCheckFunctionNorm(snes, fnorm);
225:     if (lssucceed) {
226:       if (snes->stol * xnorm > ynorm) {
227:         snes->reason = SNES_CONVERGED_SNORM_RELATIVE;
228:         return 0;
229:       }
230:       if (++snes->numFailures >= snes->maxFailures) {
231:         PetscBool ismin;
232:         snes->reason = SNES_DIVERGED_LINE_SEARCH;
233:         SNESNEWTONLSCheckLocalMin_Private(snes, snes->jacobian, F, fnorm, &ismin);
234:         if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
235:         if (snes->errorifnotconverged && snes->reason) {
236:           PetscViewer monitor;
237:           SNESLineSearchGetDefaultMonitor(linesearch, &monitor);
239:           SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_NOT_CONVERGED, "SNESSolve has not converged due %s.", SNESConvergedReasons[snes->reason]);
240:         }
241:         break;
242:       }
243:     }
244:     /* Monitor convergence */
245:     PetscObjectSAWsTakeAccess((PetscObject)snes);
246:     snes->iter  = i + 1;
247:     snes->norm  = fnorm;
248:     snes->ynorm = ynorm;
249:     snes->xnorm = xnorm;
250:     PetscObjectSAWsGrantAccess((PetscObject)snes);
251:     SNESLogConvergenceHistory(snes, snes->norm, lits);
252:     SNESMonitor(snes, snes->iter, snes->norm);
253:     /* Test for convergence */
254:     PetscUseTypeMethod(snes, converged, snes->iter, xnorm, ynorm, fnorm, &snes->reason, snes->cnvP);
255:     if (snes->reason) break;
256:   }
257:   if (i == maxits) {
258:     PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", maxits);
259:     if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
260:   }
261:   return 0;
262: }

264: /*
265:    SNESSetUp_NEWTONLS - Sets up the internal data structures for the later use
266:    of the SNESNEWTONLS nonlinear solver.

268:    Input Parameter:
269: .  snes - the SNES context
270: .  x - the solution vector

272:    Application Interface Routine: SNESSetUp()

274:  */
275: PetscErrorCode SNESSetUp_NEWTONLS(SNES snes)
276: {
277:   SNESSetUpMatrices(snes);
278:   if (snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_DEFAULT) snes->functype = SNES_FUNCTION_PRECONDITIONED;
279:   return 0;
280: }

282: PetscErrorCode SNESReset_NEWTONLS(SNES snes)
283: {
284:   return 0;
285: }

287: /*
288:    SNESDestroy_NEWTONLS - Destroys the private SNES_NEWTONLS context that was created
289:    with SNESCreate_NEWTONLS().

291:    Input Parameter:
292: .  snes - the SNES context

294:    Application Interface Routine: SNESDestroy()
295:  */
296: PetscErrorCode SNESDestroy_NEWTONLS(SNES snes)
297: {
298:   SNESReset_NEWTONLS(snes);
299:   PetscFree(snes->data);
300:   return 0;
301: }

303: /*
304:    SNESView_NEWTONLS - Prints info from the SNESNEWTONLS data structure.

306:    Input Parameters:
307: .  SNES - the SNES context
308: .  viewer - visualization context

310:    Application Interface Routine: SNESView()
311: */
312: static PetscErrorCode SNESView_NEWTONLS(SNES snes, PetscViewer viewer)
313: {
314:   PetscBool iascii;

316:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
317:   if (iascii) { }
318:   return 0;
319: }

321: /*
322:    SNESSetFromOptions_NEWTONLS - Sets various parameters for the SNESNEWTONLS method.

324:    Input Parameter:
325: .  snes - the SNES context

327:    Application Interface Routine: SNESSetFromOptions()
328: */
329: static PetscErrorCode SNESSetFromOptions_NEWTONLS(SNES snes, PetscOptionItems *PetscOptionsObject)
330: {
331:   return 0;
332: }

334: /*MC
335:       SNESNEWTONLS - Newton based nonlinear solver that uses a line search

337:    Options Database Keys:
338: +   -snes_linesearch_type <bt> - bt,basic.  Select line search type
339: .   -snes_linesearch_order <3> - 2, 3. Selects the order of the line search for bt
340: .   -snes_linesearch_norms <true> - Turns on/off computation of the norms for basic linesearch (`SNESLineSearchSetComputeNorms()`)
341: .   -snes_linesearch_alpha <alpha> - Sets alpha used in determining if reduction in function norm is sufficient
342: .   -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)
343: .   -snes_linesearch_minlambda <minlambda>  - Sets the minimum lambda the line search will tolerate
344: .   -snes_linesearch_monitor - print information about progress of line searches
345: -   -snes_linesearch_damping - damping factor used for basic line search

347:     Note:
348:     This is the default nonlinear solver in `SNES`

350:    Level: beginner

352: .seealso: `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONTR`, `SNESQN`, `SNESLineSearchSetType()`, `SNESLineSearchSetOrder()`
353:           `SNESLineSearchSetPostCheck()`, `SNESLineSearchSetPreCheck()` `SNESLineSearchSetComputeNorms()`, `SNESGetLineSearch()`
354: M*/
355: PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONLS(SNES snes)
356: {
357:   SNES_NEWTONLS *neP;
358:   SNESLineSearch linesearch;

360:   snes->ops->setup          = SNESSetUp_NEWTONLS;
361:   snes->ops->solve          = SNESSolve_NEWTONLS;
362:   snes->ops->destroy        = SNESDestroy_NEWTONLS;
363:   snes->ops->setfromoptions = SNESSetFromOptions_NEWTONLS;
364:   snes->ops->view           = SNESView_NEWTONLS;
365:   snes->ops->reset          = SNESReset_NEWTONLS;

367:   snes->npcside = PC_RIGHT;
368:   snes->usesksp = PETSC_TRUE;
369:   snes->usesnpc = PETSC_TRUE;

371:   SNESGetLineSearch(snes, &linesearch);
372:   if (!((PetscObject)linesearch)->type_name) SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);

374:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

376:   PetscNew(&neP);
377:   snes->data = (void *)neP;
378:   return 0;
379: }