Actual source code: ls.c

petsc-3.7.4 2016-10-02
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: */
 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;

151:   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);

153:   snes->numFailures            = 0;
154:   snes->numLinearSolveFailures = 0;
155:   snes->reason                 = SNES_CONVERGED_ITERATING;

157:   maxits = snes->max_its;               /* maximum number of iterations */
158:   X      = snes->vec_sol;               /* solution vector */
159:   F      = snes->vec_func;              /* residual vector */
160:   Y      = snes->vec_sol_update;        /* newton step */

162:   PetscObjectSAWsTakeAccess((PetscObject)snes);
163:   snes->iter = 0;
164:   snes->norm = 0.0;
165:   PetscObjectSAWsGrantAccess((PetscObject)snes);
166:   SNESGetLineSearch(snes, &linesearch);

168:   /* compute the preconditioned function first in the case of left preconditioning with preconditioned function */
169:   if (snes->pc && snes->pcside == PC_LEFT && snes->functype == SNES_FUNCTION_PRECONDITIONED) {
170:     SNESApplyNPC(snes,X,NULL,F);
171:     SNESGetConvergedReason(snes->pc,&reason);
172:     if (reason < 0  && reason != SNES_DIVERGED_MAX_IT) {
173:       snes->reason = SNES_DIVERGED_INNER;
174:       return(0);
175:     }

177:     VecNormBegin(F,NORM_2,&fnorm);
178:     VecNormEnd(F,NORM_2,&fnorm);
179:   } else {
180:     if (!snes->vec_func_init_set) {
181:       SNESComputeFunction(snes,X,F);
182:     } else snes->vec_func_init_set = PETSC_FALSE;
183:   }

185:   VecNorm(F,NORM_2,&fnorm);        /* fnorm <- ||F||  */
186:   SNESCheckFunctionNorm(snes,fnorm);
187:   PetscObjectSAWsTakeAccess((PetscObject)snes);
188:   snes->norm = fnorm;
189:   PetscObjectSAWsGrantAccess((PetscObject)snes);
190:   SNESLogConvergenceHistory(snes,fnorm,0);
191:   SNESMonitor(snes,0,fnorm);

193:   /* test convergence */
194:   (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
195:   if (snes->reason) return(0);

197:   for (i=0; i<maxits; i++) {

199:     /* Call general purpose update function */
200:     if (snes->ops->update) {
201:       (*snes->ops->update)(snes, snes->iter);
202:     }

204:     /* apply the nonlinear preconditioner */
205:     if (snes->pc) {
206:       if (snes->pcside == PC_RIGHT) {
207:         SNESSetInitialFunction(snes->pc, F);
208:         PetscLogEventBegin(SNES_NPCSolve,snes->pc,X,snes->vec_rhs,0);
209:         SNESSolve(snes->pc, snes->vec_rhs, X);
210:         PetscLogEventEnd(SNES_NPCSolve,snes->pc,X,snes->vec_rhs,0);
211:         SNESGetConvergedReason(snes->pc,&reason);
212:         if (reason < 0  && reason != SNES_DIVERGED_MAX_IT) {
213:           snes->reason = SNES_DIVERGED_INNER;
214:           return(0);
215:         }
216:         SNESGetNPCFunction(snes,F,&fnorm);
217:       } else if (snes->pcside == PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
218:         SNESApplyNPC(snes,X,F,F);
219:         SNESGetConvergedReason(snes->pc,&reason);
220:         if (reason < 0  && reason != SNES_DIVERGED_MAX_IT) {
221:           snes->reason = SNES_DIVERGED_INNER;
222:           return(0);
223:         }
224:       }
225:     }

227:     /* Solve J Y = F, where J is Jacobian matrix */
228:     SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
229:     KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);
230:     KSPSolve(snes->ksp,F,Y);
231:     SNESCheckKSPSolve(snes);
232:     KSPGetIterationNumber(snes->ksp,&lits);
233:     snes->linear_its += lits;
234:     PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);

236:     if (PetscLogPrintInfo) {
237:       SNESNEWTONLSCheckResidual_Private(snes,snes->jacobian,F,Y);
238:     }

240:     /* Compute a (scaled) negative update in the line search routine:
241:          X <- X - lambda*Y
242:        and evaluate F = function(X) (depends on the line search).
243:     */
244:     gnorm = fnorm;
245:     SNESLineSearchApply(linesearch, X, F, &fnorm, Y);
246:     SNESLineSearchGetReason(linesearch, &lssucceed);
247:     SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm);
248:     PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",(double)gnorm,(double)fnorm,(double)ynorm,(int)lssucceed);
249:     if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
250:     SNESCheckFunctionNorm(snes,fnorm);
251:     if (lssucceed) {
252:       if (snes->stol*xnorm > ynorm) {
253:         snes->reason = SNES_CONVERGED_SNORM_RELATIVE;
254:         return(0);
255:       }
256:       if (++snes->numFailures >= snes->maxFailures) {
257:         PetscBool ismin;
258:         snes->reason = SNES_DIVERGED_LINE_SEARCH;
259:         SNESNEWTONLSCheckLocalMin_Private(snes,snes->jacobian,F,fnorm,&ismin);
260:         if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
261:         break;
262:       }
263:     }
264:     /* Monitor convergence */
265:     PetscObjectSAWsTakeAccess((PetscObject)snes);
266:     snes->iter = i+1;
267:     snes->norm = fnorm;
268:     PetscObjectSAWsGrantAccess((PetscObject)snes);
269:     SNESLogConvergenceHistory(snes,snes->norm,lits);
270:     SNESMonitor(snes,snes->iter,snes->norm);
271:     /* Test for convergence */
272:     (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);
273:     if (snes->reason) break;
274:   }
275:   if (i == maxits) {
276:     PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);
277:     if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
278:   }
279:   return(0);
280: }
281: /* -------------------------------------------------------------------------- */
282: /*
283:    SNESSetUp_NEWTONLS - Sets up the internal data structures for the later use
284:    of the SNESNEWTONLS nonlinear solver.

286:    Input Parameter:
287: .  snes - the SNES context
288: .  x - the solution vector

290:    Application Interface Routine: SNESSetUp()

292:    Notes:
293:    For basic use of the SNES solvers, the user need not explicitly call
294:    SNESSetUp(), since these actions will automatically occur during
295:    the call to SNESSolve().
296:  */
299: PetscErrorCode SNESSetUp_NEWTONLS(SNES snes)
300: {

304:   SNESSetUpMatrices(snes);
305:   if (snes->pcside == PC_LEFT && snes->functype == SNES_FUNCTION_DEFAULT) snes->functype = SNES_FUNCTION_PRECONDITIONED;
306:   return(0);
307: }
308: /* -------------------------------------------------------------------------- */

312: PetscErrorCode SNESReset_NEWTONLS(SNES snes)
313: {
315:   return(0);
316: }

318: /*
319:    SNESDestroy_NEWTONLS - Destroys the private SNES_NEWTONLS context that was created
320:    with SNESCreate_NEWTONLS().

322:    Input Parameter:
323: .  snes - the SNES context

325:    Application Interface Routine: SNESDestroy()
326:  */
329: PetscErrorCode SNESDestroy_NEWTONLS(SNES snes)
330: {

334:   SNESReset_NEWTONLS(snes);
335:   PetscFree(snes->data);
336:   return(0);
337: }
338: /* -------------------------------------------------------------------------- */

340: /*
341:    SNESView_NEWTONLS - Prints info from the SNESNEWTONLS data structure.

343:    Input Parameters:
344: .  SNES - the SNES context
345: .  viewer - visualization context

347:    Application Interface Routine: SNESView()
348: */
351: static PetscErrorCode SNESView_NEWTONLS(SNES snes,PetscViewer viewer)
352: {
354:   PetscBool      iascii;

357:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
358:   if (iascii) {
359:   }
360:   return(0);
361: }

363: /* -------------------------------------------------------------------------- */
364: /*
365:    SNESSetFromOptions_NEWTONLS - Sets various parameters for the SNESNEWTONLS method.

367:    Input Parameter:
368: .  snes - the SNES context

370:    Application Interface Routine: SNESSetFromOptions()
371: */
374: static PetscErrorCode SNESSetFromOptions_NEWTONLS(PetscOptionItems *PetscOptionsObject,SNES snes)
375: {
377:   SNESLineSearch linesearch;

380:   if (!snes->linesearch) {
381:     SNESGetLineSearch(snes, &linesearch);
382:     SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);
383:   }
384:   return(0);
385: }

387: /* -------------------------------------------------------------------------- */
388: /*MC
389:       SNESNEWTONLS - Newton based nonlinear solver that uses a line search

391:    Options Database:
392: +   -snes_linesearch_type <bt> - bt,basic.  Select line search type
393: .   -snes_linesearch_order <3> - 2, 3. Selects the order of the line search for bt
394: .   -snes_linesearch_norms <true> - Turns on/off computation of the norms for basic linesearch
395: .   -snes_linesearch_alpha <alpha> - Sets alpha used in determining if reduction in function norm is sufficient
396: .   -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)
397: .   -snes_linesearch_minlambda <minlambda>  - Sets the minimum lambda the line search will tolerate
398: .   -snes_linesearch_monitor - print information about progress of line searches
399: -   -snes_linesearch_damping - damping factor used for basic line search

401:     Notes: This is the default nonlinear solver in SNES

403:    Level: beginner

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

408: M*/
411: PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONLS(SNES snes)
412: {
414:   SNES_NEWTONLS  *neP;

417:   snes->ops->setup          = SNESSetUp_NEWTONLS;
418:   snes->ops->solve          = SNESSolve_NEWTONLS;
419:   snes->ops->destroy        = SNESDestroy_NEWTONLS;
420:   snes->ops->setfromoptions = SNESSetFromOptions_NEWTONLS;
421:   snes->ops->view           = SNESView_NEWTONLS;
422:   snes->ops->reset          = SNESReset_NEWTONLS;

424:   snes->pcside  = PC_RIGHT;
425:   snes->usesksp = PETSC_TRUE;
426:   snes->usespc  = PETSC_TRUE;
427:   PetscNewLog(snes,&neP);
428:   snes->data    = (void*)neP;
429:   return(0);
430: }