Actual source code: ls.c
1: #define PETSCSNES_DLL
3: #include src/snes/impls/ls/ls.h
5: /*
6: Checks if J^T F = 0 which implies we've found a local minimum of the function,
7: but not a zero. In the case when one cannot compute J^T F we use the fact that
8: 0 = (J^T F)^T W = F^T J W iff W not in the null space of J. Thanks for Jorge More
9: for this trick.
10: */
13: PetscErrorCode SNESLSCheckLocalMin_Private(Mat A,Vec F,Vec W,PetscReal fnorm,PetscTruth *ismin)
14: {
15: PetscReal a1;
17: PetscTruth hastranspose;
20: *ismin = PETSC_FALSE;
21: MatHasOperation(A,MATOP_MULT_TRANSPOSE,&hastranspose);
22: if (hastranspose) {
23: /* Compute || J^T F|| */
24: MatMultTranspose(A,F,W);
25: VecNorm(W,NORM_2,&a1);
26: PetscInfo1(0,"|| J^T F|| %G near zero implies found a local minimum\n",a1/fnorm);
27: if (a1/fnorm < 1.e-4) *ismin = PETSC_TRUE;
28: } else {
29: Vec work;
30: PetscScalar result;
31: PetscReal wnorm;
33: VecSetRandom(W,PETSC_NULL);
34: VecNorm(W,NORM_2,&wnorm);
35: VecDuplicate(W,&work);
36: MatMult(A,W,work);
37: VecDot(F,work,&result);
38: VecDestroy(work);
39: a1 = PetscAbsScalar(result)/(fnorm*wnorm);
40: PetscInfo1(0,"(F^T J random)/(|| F ||*||random|| %G near zero implies found a local minimum\n",a1);
41: if (a1 < 1.e-4) *ismin = PETSC_TRUE;
42: }
43: return(0);
44: }
46: /*
47: Checks if J^T(F - J*X) = 0
48: */
51: PetscErrorCode SNESLSCheckResidual_Private(Mat A,Vec F,Vec X,Vec W1,Vec W2)
52: {
53: PetscReal a1,a2;
55: PetscTruth hastranspose;
58: MatHasOperation(A,MATOP_MULT_TRANSPOSE,&hastranspose);
59: if (hastranspose) {
60: MatMult(A,X,W1);
61: VecAXPY(W1,-1.0,F);
63: /* Compute || J^T W|| */
64: MatMultTranspose(A,W1,W2);
65: VecNorm(W1,NORM_2,&a1);
66: VecNorm(W2,NORM_2,&a2);
67: if (a1) {
68: PetscInfo1(0,"||J^T(F-Ax)||/||F-AX|| %G near zero implies inconsistent rhs\n",a2/a1);
69: }
70: }
71: return(0);
72: }
74: /* --------------------------------------------------------------------
76: This file implements a truncated Newton method with a line search,
77: for solving a system of nonlinear equations, using the KSP, Vec,
78: and Mat interfaces for linear solvers, vectors, and matrices,
79: respectively.
81: The following basic routines are required for each nonlinear solver:
82: SNESCreate_XXX() - Creates a nonlinear solver context
83: SNESSetFromOptions_XXX() - Sets runtime options
84: SNESSolve_XXX() - Solves the nonlinear system
85: SNESDestroy_XXX() - Destroys the nonlinear solver context
86: The suffix "_XXX" denotes a particular implementation, in this case
87: we use _LS (e.g., SNESCreate_LS, SNESSolve_LS) for solving
88: systems of nonlinear equations with a line search (LS) method.
89: These routines are actually called via the common user interface
90: routines SNESCreate(), SNESSetFromOptions(), SNESSolve(), and
91: SNESDestroy(), so the application code interface remains identical
92: for all nonlinear solvers.
94: Another key routine is:
95: SNESSetUp_XXX() - Prepares for the use of a nonlinear solver
96: by setting data structures and options. The interface routine SNESSetUp()
97: is not usually called directly by the user, but instead is called by
98: SNESSolve() if necessary.
100: Additional basic routines are:
101: SNESView_XXX() - Prints details of runtime options that
102: have actually been used.
103: These are called by application codes via the interface routines
104: SNESView().
106: The various types of solvers (preconditioners, Krylov subspace methods,
107: nonlinear solvers, timesteppers) are all organized similarly, so the
108: above description applies to these categories also.
110: -------------------------------------------------------------------- */
111: /*
112: SNESSolve_LS - Solves a nonlinear system with a truncated Newton
113: method with a line search.
115: Input Parameters:
116: . snes - the SNES context
118: Output Parameter:
119: . outits - number of iterations until termination
121: Application Interface Routine: SNESSolve()
123: Notes:
124: This implements essentially a truncated Newton method with a
125: line search. By default a cubic backtracking line search
126: is employed, as described in the text "Numerical Methods for
127: Unconstrained Optimization and Nonlinear Equations" by Dennis
128: and Schnabel.
129: */
132: PetscErrorCode SNESSolve_LS(SNES snes)
133: {
134: SNES_LS *neP = (SNES_LS*)snes->data;
135: PetscErrorCode ierr;
136: PetscInt maxits,i,lits;
137: PetscTruth lssucceed;
138: MatStructure flg = DIFFERENT_NONZERO_PATTERN;
139: PetscReal fnorm,gnorm,xnorm,ynorm;
140: Vec Y,X,F,G,W,TMP;
141: KSPConvergedReason kspreason;
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->work[0]; /* work vectors */
152: G = snes->work[1];
153: W = snes->work[2];
155: PetscObjectTakeAccess(snes);
156: snes->iter = 0;
157: PetscObjectGrantAccess(snes);
158: SNESComputeFunction(snes,X,F);
159: VecNorm(F,NORM_2,&fnorm); /* fnorm <- ||F|| */
160: if (fnorm != fnorm) SETERRQ(PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
161: PetscObjectTakeAccess(snes);
162: snes->norm = fnorm;
163: PetscObjectGrantAccess(snes);
164: SNESLogConvHistory(snes,fnorm,0);
165: SNESMonitor(snes,0,fnorm);
166: (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
167: if (snes->reason) return(0);
169: for (i=0; i<maxits; i++) {
171: /* Call general purpose update function */
172: if (snes->ops->update) {
173: (*snes->ops->update)(snes, snes->iter);
174: }
176: /* Solve J Y = F, where J is Jacobian matrix */
177: SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);
178: KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,flg);
179: KSPSolve(snes->ksp,F,Y);
180: KSPGetConvergedReason(snes->ksp,&kspreason);
181: if (kspreason < 0) {
182: if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
183: snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
184: return(0);
185: }
186: }
187: KSPGetIterationNumber(snes->ksp,&lits);
189: if (neP->precheckstep) {
190: PetscTruth changed_y = PETSC_FALSE;
191: (*neP->precheckstep)(snes,X,Y,neP->precheck,&changed_y);
192: }
194: if (PetscLogPrintInfo){
195: SNESLSCheckResidual_Private(snes->jacobian,F,Y,G,W);
196: }
198: /* should check what happened to the linear solve? */
199: snes->linear_its += lits;
200: PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);
202: /* Compute a (scaled) negative update in the line search routine:
203: Y <- X - lambda*Y
204: and evaluate G = function(Y) (depends on the line search).
205: */
206: VecCopy(Y,snes->vec_sol_update_always);
207: (*neP->LineSearch)(snes,neP->lsP,X,F,G,Y,W,fnorm,&ynorm,&gnorm,&lssucceed);
208: PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",fnorm,gnorm,ynorm,(int)lssucceed);
209: TMP = F; F = G; snes->vec_func_always = F; G = TMP;
210: TMP = X; X = W; snes->vec_sol_always = X; W = TMP;
211: fnorm = gnorm;
212: if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
214: PetscObjectTakeAccess(snes);
215: snes->iter = i+1;
216: snes->norm = fnorm;
217: PetscObjectGrantAccess(snes);
218: SNESLogConvHistory(snes,fnorm,lits);
219: SNESMonitor(snes,i+1,fnorm);
221: if (!lssucceed) {
222: PetscTruth ismin;
224: if (++snes->numFailures >= snes->maxFailures) {
225: snes->reason = SNES_DIVERGED_LS_FAILURE;
226: SNESLSCheckLocalMin_Private(snes->jacobian,F,W,fnorm,&ismin);
227: if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
228: break;
229: }
230: }
232: /* Test for convergence */
233: if (snes->ops->converged) {
234: VecNorm(X,NORM_2,&xnorm); /* xnorm = || X || */
235: (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);
236: if (snes->reason) {
237: break;
238: }
239: }
240: }
241: if (X != snes->vec_sol) {
242: VecCopy(X,snes->vec_sol);
243: }
244: if (F != snes->vec_func) {
245: VecCopy(F,snes->vec_func);
246: }
247: snes->vec_sol_always = snes->vec_sol;
248: snes->vec_func_always = snes->vec_func;
249: if (i == maxits) {
250: PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);
251: snes->reason = SNES_DIVERGED_MAX_IT;
252: }
253: return(0);
254: }
255: /* -------------------------------------------------------------------------- */
256: /*
257: SNESSetUp_LS - Sets up the internal data structures for the later use
258: of the SNESLS nonlinear solver.
260: Input Parameter:
261: . snes - the SNES context
262: . x - the solution vector
264: Application Interface Routine: SNESSetUp()
266: Notes:
267: For basic use of the SNES solvers, the user need not explicitly call
268: SNESSetUp(), since these actions will automatically occur during
269: the call to SNESSolve().
270: */
273: PetscErrorCode SNESSetUp_LS(SNES snes)
274: {
278: if (!snes->work) {
279: snes->nwork = 4;
280: VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work);
281: PetscLogObjectParents(snes,snes->nwork,snes->work);
282: snes->vec_sol_update_always = snes->work[3];
283: }
284: return(0);
285: }
286: /* -------------------------------------------------------------------------- */
287: /*
288: SNESDestroy_LS - Destroys the private SNES_LS context that was created
289: with SNESCreate_LS().
291: Input Parameter:
292: . snes - the SNES context
294: Application Interface Routine: SNESDestroy()
295: */
298: PetscErrorCode SNESDestroy_LS(SNES snes)
299: {
303: if (snes->nwork) {
304: VecDestroyVecs(snes->work,snes->nwork);
305: }
306: PetscFree(snes->data);
307: return(0);
308: }
309: /* -------------------------------------------------------------------------- */
313: /*@C
314: SNESLineSearchNo - This routine is not a line search at all;
315: it simply uses the full Newton step. Thus, this routine is intended
316: to serve as a template and is not recommended for general use.
318: Collective on SNES and Vec
320: Input Parameters:
321: + snes - nonlinear context
322: . lsctx - optional context for line search (not used here)
323: . x - current iterate
324: . f - residual evaluated at x
325: . y - search direction
326: . w - work vector
327: - fnorm - 2-norm of f
329: Output Parameters:
330: + g - residual evaluated at new iterate y
331: . w - new iterate
332: . gnorm - 2-norm of g
333: . ynorm - 2-norm of search length
334: - flag - PETSC_TRUE on success, PETSC_FALSE on failure
336: Options Database Key:
337: . -snes_ls basic - Activates SNESLineSearchNo()
339: Level: advanced
341: .keywords: SNES, nonlinear, line search, cubic
343: .seealso: SNESLineSearchCubic(), SNESLineSearchQuadratic(),
344: SNESLineSearchSet(), SNESLineSearchNoNorms()
345: @*/
346: PetscErrorCode SNESLineSearchNo(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal *ynorm,PetscReal *gnorm,PetscTruth *flag)
347: {
349: SNES_LS *neP = (SNES_LS*)snes->data;
350: PetscTruth changed_w = PETSC_FALSE,changed_y = PETSC_FALSE;
353: *flag = PETSC_TRUE;
355: VecNorm(y,NORM_2,ynorm); /* ynorm = || y || */
356: VecWAXPY(w,-1.0,y,x); /* w <- x - y */
357: if (neP->postcheckstep) {
358: (*neP->postcheckstep)(snes,x,y,w,neP->postcheck,&changed_y,&changed_w);
359: }
360: if (changed_y) {
361: VecWAXPY(w,-1.0,y,x); /* w <- x - y */
362: }
363: SNESComputeFunction(snes,w,g);
364: if (PetscExceptionValue(ierr)) {
366: }
367:
369: VecNorm(g,NORM_2,gnorm); /* gnorm = || g || */
370: if (*gnorm != *gnorm) SETERRQ(PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
372: return(0);
373: }
374: /* -------------------------------------------------------------------------- */
378: /*@C
379: SNESLineSearchNoNorms - This routine is not a line search at
380: all; it simply uses the full Newton step. This version does not
381: even compute the norm of the function or search direction; this
382: is intended only when you know the full step is fine and are
383: not checking for convergence of the nonlinear iteration (for
384: example, you are running always for a fixed number of Newton steps).
386: Collective on SNES and Vec
388: Input Parameters:
389: + snes - nonlinear context
390: . lsctx - optional context for line search (not used here)
391: . x - current iterate
392: . f - residual evaluated at x
393: . y - search direction
394: . w - work vector
395: - fnorm - 2-norm of f
397: Output Parameters:
398: + g - residual evaluated at new iterate y
399: . w - new iterate
400: . gnorm - not changed
401: . ynorm - not changed
402: - flag - set to PETSC_TRUE indicating a successful line search
404: Options Database Key:
405: . -snes_ls basicnonorms - Activates SNESLineSearchNoNorms()
407: Notes:
408: SNESLineSearchNoNorms() must be used in conjunction with
409: either the options
410: $ -snes_no_convergence_test -snes_max_it <its>
411: or alternatively a user-defined custom test set via
412: SNESSetConvergenceTest(); or a -snes_max_it of 1,
413: otherwise, the SNES solver will generate an error.
415: During the final iteration this will not evaluate the function at
416: the solution point. This is to save a function evaluation while
417: using pseudo-timestepping.
419: The residual norms printed by monitoring routines such as
420: SNESMonitorDefault() (as activated via -snes_monitor) will not be
421: correct, since they are not computed.
423: Level: advanced
425: .keywords: SNES, nonlinear, line search, cubic
427: .seealso: SNESLineSearchCubic(), SNESLineSearchQuadratic(),
428: SNESLineSearchSet(), SNESLineSearchNo()
429: @*/
430: PetscErrorCode SNESLineSearchNoNorms(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal *ynorm,PetscReal *gnorm,PetscTruth *flag)
431: {
433: SNES_LS *neP = (SNES_LS*)snes->data;
434: PetscTruth changed_w = PETSC_FALSE,changed_y = PETSC_FALSE;
437: *flag = PETSC_TRUE;
439: VecWAXPY(w,-1.0,y,x); /* w <- x - y */
440: if (neP->postcheckstep) {
441: (*neP->postcheckstep)(snes,x,y,w,neP->postcheck,&changed_y,&changed_w);
442: }
443: if (changed_y) {
444: VecWAXPY(w,-1.0,y,x); /* w <- x - y */
445: }
446:
447: /* don't evaluate function the last time through */
448: if (snes->iter < snes->max_its-1) {
449: SNESComputeFunction(snes,w,g);
450: if (PetscExceptionValue(ierr)) {
452: }
453:
454: }
456: return(0);
457: }
458: /* -------------------------------------------------------------------------- */
461: /*@C
462: SNESLineSearchCubic - Performs a cubic line search (default line search method).
464: Collective on SNES
466: Input Parameters:
467: + snes - nonlinear context
468: . lsctx - optional context for line search (not used here)
469: . x - current iterate
470: . f - residual evaluated at x
471: . y - search direction
472: . w - work vector
473: - fnorm - 2-norm of f
475: Output Parameters:
476: + g - residual evaluated at new iterate y
477: . w - new iterate
478: . gnorm - 2-norm of g
479: . ynorm - 2-norm of search length
480: - flag - PETSC_TRUE if line search succeeds; PETSC_FALSE on failure.
482: Options Database Key:
483: $ -snes_ls cubic - Activates SNESLineSearchCubic()
485: Notes:
486: This line search is taken from "Numerical Methods for Unconstrained
487: Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325.
489: Level: advanced
491: .keywords: SNES, nonlinear, line search, cubic
493: .seealso: SNESLineSearchQuadratic(), SNESLineSearchNo(), SNESLineSearchSet(), SNESLineSearchNoNorms()
494: @*/
495: PetscErrorCode SNESLineSearchCubic(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal *ynorm,PetscReal *gnorm,PetscTruth *flag)
496: {
497: /*
498: Note that for line search purposes we work with with the related
499: minimization problem:
500: min z(x): R^n -> R,
501: where z(x) = .5 * fnorm*fnorm, and fnorm = || f ||_2.
502: */
503:
504: PetscReal steptol,initslope,lambdaprev,gnormprev,a,b,d,t1,t2,rellength;
505: PetscReal maxstep,minlambda,alpha,lambda,lambdatemp;
506: #if defined(PETSC_USE_COMPLEX)
507: PetscScalar cinitslope,clambda;
508: #endif
510: PetscInt count;
511: SNES_LS *neP = (SNES_LS*)snes->data;
512: PetscScalar scale;
513: PetscTruth changed_w = PETSC_FALSE,changed_y = PETSC_FALSE;
517: *flag = PETSC_TRUE;
518: alpha = neP->alpha;
519: maxstep = neP->maxstep;
520: steptol = neP->steptol;
522: VecNorm(y,NORM_2,ynorm);
523: if (!*ynorm) {
524: PetscInfo(snes,"Search direction and size is 0\n");
525: *gnorm = fnorm;
526: VecCopy(x,w);
527: VecCopy(f,g);
528: *flag = PETSC_FALSE;
529: goto theend1;
530: }
531: if (*ynorm > maxstep) { /* Step too big, so scale back */
532: scale = maxstep/(*ynorm);
533: PetscInfo2(snes,"Scaling step by %G old ynorm %G\n",PetscRealPart(scale),*ynorm);
534: VecScale(y,scale);
535: *ynorm = maxstep;
536: }
537: VecMaxPointwiseDivide(y,x,&rellength);
538: minlambda = steptol/rellength;
539: MatMult(snes->jacobian,y,w);
540: #if defined(PETSC_USE_COMPLEX)
541: VecDot(f,w,&cinitslope);
542: initslope = PetscRealPart(cinitslope);
543: #else
544: VecDot(f,w,&initslope);
545: #endif
546: if (initslope > 0.0) initslope = -initslope;
547: if (initslope == 0.0) initslope = -1.0;
549: VecWAXPY(w,-1.0,y,x);
550: if (snes->nfuncs >= snes->max_funcs) {
551: PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");
552: *flag = PETSC_FALSE;
553: snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
554: goto theend1;
555: }
556: SNESComputeFunction(snes,w,g);
557: if (PetscExceptionValue(ierr)) {
559: }
560:
561: VecNorm(g,NORM_2,gnorm);
562: if (*gnorm != *gnorm) SETERRQ(PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
563: if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + alpha*initslope) { /* Sufficient reduction */
564: PetscInfo(snes,"Using full step\n");
565: goto theend1;
566: }
568: /* Fit points with quadratic */
569: lambda = 1.0;
570: lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope);
571: lambdaprev = lambda;
572: gnormprev = *gnorm;
573: if (lambdatemp > .5*lambda) lambdatemp = .5*lambda;
574: if (lambdatemp <= .1*lambda) lambda = .1*lambda;
575: else lambda = lambdatemp;
577: #if defined(PETSC_USE_COMPLEX)
578: clambda = lambda; VecWAXPY(w,-clambda,y,x);
579: #else
580: VecWAXPY(w,-lambda,y,x);
581: #endif
582: if (snes->nfuncs >= snes->max_funcs) {
583: PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);
584: *flag = PETSC_FALSE;
585: snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
586: goto theend1;
587: }
588: SNESComputeFunction(snes,w,g);
589: if (PetscExceptionValue(ierr)) {
591: }
592:
593: VecNorm(g,NORM_2,gnorm);
594: if (*gnorm != *gnorm) SETERRQ(PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
595: if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*alpha*initslope) { /* sufficient reduction */
596: PetscInfo1(snes,"Quadratically determined step, lambda=%18.16e\n",lambda);
597: goto theend1;
598: }
600: /* Fit points with cubic */
601: count = 1;
602: while (count < 10000) {
603: if (lambda <= minlambda) { /* bad luck; use full step */
604: PetscInfo1(snes,"Unable to find good step length! %D \n",count);
605: PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",fnorm,*gnorm,*ynorm,lambda,initslope);
606: *flag = PETSC_FALSE;
607: break;
608: }
609: t1 = .5*((*gnorm)*(*gnorm) - fnorm*fnorm) - lambda*initslope;
610: t2 = .5*(gnormprev*gnormprev - fnorm*fnorm) - lambdaprev*initslope;
611: a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
612: b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
613: d = b*b - 3*a*initslope;
614: if (d < 0.0) d = 0.0;
615: if (a == 0.0) {
616: lambdatemp = -initslope/(2.0*b);
617: } else {
618: lambdatemp = (-b + sqrt(d))/(3.0*a);
619: }
620: lambdaprev = lambda;
621: gnormprev = *gnorm;
622: if (lambdatemp > .5*lambda) lambdatemp = .5*lambda;
623: if (lambdatemp <= .1*lambda) lambda = .1*lambda;
624: else lambda = lambdatemp;
625: #if defined(PETSC_USE_COMPLEX)
626: clambda = lambda;
627: VecWAXPY(w,-clambda,y,x);
628: #else
629: VecWAXPY(w,-lambda,y,x);
630: #endif
631: if (snes->nfuncs >= snes->max_funcs) {
632: PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);
633: PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",fnorm,*gnorm,*ynorm,lambda,initslope);
634: *flag = PETSC_FALSE;
635: snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
636: break;
637: }
638: SNESComputeFunction(snes,w,g);
639: if (PetscExceptionValue(ierr)) {
641: }
642:
643: VecNorm(g,NORM_2,gnorm);
644: if (*gnorm != *gnorm) SETERRQ(PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
645: if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*alpha*initslope) { /* is reduction enough? */
646: PetscInfo1(snes,"Cubically determined step, lambda=%18.16e\n",lambda);
647: break;
648: } else {
649: PetscInfo1(snes,"Cubic step no good, shrinking lambda, lambda=%18.16e\n",lambda);
650: }
651: count++;
652: }
653: if (count >= 10000) {
654: SETERRQ(PETSC_ERR_LIB, "Lambda was decreased more than 10,000 times, so something is probably wrong with the function evaluation");
655: }
656: theend1:
657: /* Optional user-defined check for line search step validity */
658: if (neP->postcheckstep && *flag) {
659: (*neP->postcheckstep)(snes,x,y,w,neP->postcheck,&changed_y,&changed_w);
660: if (changed_y) {
661: VecWAXPY(w,-1.0,y,x);
662: }
663: if (changed_y || changed_w) { /* recompute the function if the step has changed */
664: SNESComputeFunction(snes,w,g);
665: if (PetscExceptionValue(ierr)) {
667: }
668:
669: VecNormBegin(g,NORM_2,gnorm);
670: if (*gnorm != *gnorm) SETERRQ(PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
671: VecNormBegin(w,NORM_2,ynorm);
672: VecNormEnd(g,NORM_2,gnorm);
673: VecNormEnd(w,NORM_2,ynorm);
674: }
675: }
677: return(0);
678: }
679: /* -------------------------------------------------------------------------- */
682: /*@C
683: SNESLineSearchQuadratic - Performs a quadratic line search.
685: Collective on SNES and Vec
687: Input Parameters:
688: + snes - the SNES context
689: . lsctx - optional context for line search (not used here)
690: . x - current iterate
691: . f - residual evaluated at x
692: . y - search direction
693: . w - work vector
694: - fnorm - 2-norm of f
696: Output Parameters:
697: + g - residual evaluated at new iterate w
698: . w - new iterate (x + alpha*y)
699: . gnorm - 2-norm of g
700: . ynorm - 2-norm of search length
701: - flag - PETSC_TRUE if line search succeeds; PETSC_FALSE on failure.
703: Options Database Key:
704: . -snes_ls quadratic - Activates SNESLineSearchQuadratic()
706: Notes:
707: Use SNESLineSearchSet() to set this routine within the SNESLS method.
709: Level: advanced
711: .keywords: SNES, nonlinear, quadratic, line search
713: .seealso: SNESLineSearchCubic(), SNESLineSearchNo(), SNESLineSearchSet(), SNESLineSearchNoNorms()
714: @*/
715: PetscErrorCode SNESLineSearchQuadratic(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal *ynorm,PetscReal *gnorm,PetscTruth *flag)
716: {
717: /*
718: Note that for line search purposes we work with with the related
719: minimization problem:
720: min z(x): R^n -> R,
721: where z(x) = .5 * fnorm*fnorm,and fnorm = || f ||_2.
722: */
723: PetscReal steptol,initslope,maxstep,minlambda,alpha,lambda,lambdatemp,rellength;
724: #if defined(PETSC_USE_COMPLEX)
725: PetscScalar cinitslope,clambda;
726: #endif
728: PetscInt count;
729: SNES_LS *neP = (SNES_LS*)snes->data;
730: PetscScalar scale;
731: PetscTruth changed_w = PETSC_FALSE,changed_y = PETSC_FALSE;
735: *flag = PETSC_TRUE;
736: alpha = neP->alpha;
737: maxstep = neP->maxstep;
738: steptol = neP->steptol;
740: VecNorm(y,NORM_2,ynorm);
741: if (*ynorm == 0.0) {
742: PetscInfo(snes,"Search direction and size is 0\n");
743: *gnorm = fnorm;
744: VecCopy(x,w);
745: VecCopy(f,g);
746: *flag = PETSC_FALSE;
747: goto theend2;
748: }
749: if (*ynorm > maxstep) { /* Step too big, so scale back */
750: scale = maxstep/(*ynorm);
751: VecScale(y,scale);
752: *ynorm = maxstep;
753: }
754: VecMaxPointwiseDivide(y,x,&rellength);
755: minlambda = steptol/rellength;
756: MatMult(snes->jacobian,y,w);
757: #if defined(PETSC_USE_COMPLEX)
758: VecDot(f,w,&cinitslope);
759: initslope = PetscRealPart(cinitslope);
760: #else
761: VecDot(f,w,&initslope);
762: #endif
763: if (initslope > 0.0) initslope = -initslope;
764: if (initslope == 0.0) initslope = -1.0;
766: VecWAXPY(w,-1.0,y,x);
767: if (snes->nfuncs >= snes->max_funcs) {
768: PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");
769: *flag = PETSC_FALSE;
770: snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
771: goto theend2;
772: }
773: SNESComputeFunction(snes,w,g);
774: if (PetscExceptionValue(ierr)) {
776: }
777:
778: VecNorm(g,NORM_2,gnorm);
779: if (*gnorm != *gnorm) SETERRQ(PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
780: if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + alpha*initslope) { /* Sufficient reduction */
781: PetscInfo(snes,"Using full step\n");
782: goto theend2;
783: }
785: /* Fit points with quadratic */
786: lambda = 1.0;
787: count = 1;
788: while (PETSC_TRUE) {
789: if (lambda <= minlambda) { /* bad luck; use full step */
790: PetscInfo1(snes,"Unable to find good step length! %D \n",count);
791: PetscInfo5(snes,"fnorm=%G, gnorm=%G, ynorm=%G, lambda=%G, initial slope=%G\n",fnorm,*gnorm,*ynorm,lambda,initslope);
792: VecCopy(x,w);
793: *flag = PETSC_FALSE;
794: break;
795: }
796: lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope);
797: if (lambdatemp > .5*lambda) lambdatemp = .5*lambda;
798: if (lambdatemp <= .1*lambda) lambda = .1*lambda;
799: else lambda = lambdatemp;
800:
801: #if defined(PETSC_USE_COMPLEX)
802: clambda = lambda; VecWAXPY(w,-clambda,y,x);
803: #else
804: VecWAXPY(w,-lambda,y,x);
805: #endif
806: if (snes->nfuncs >= snes->max_funcs) {
807: PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);
808: PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",fnorm,*gnorm,*ynorm,lambda,initslope);
809: *flag = PETSC_FALSE;
810: snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
811: break;
812: }
813: SNESComputeFunction(snes,w,g);
814: if (PetscExceptionValue(ierr)) {
816: }
817:
818: VecNorm(g,NORM_2,gnorm);
819: if (*gnorm != *gnorm) SETERRQ(PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
820: if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*alpha*initslope) { /* sufficient reduction */
821: PetscInfo1(snes,"Quadratically determined step, lambda=%G\n",lambda);
822: break;
823: }
824: count++;
825: }
826: theend2:
827: /* Optional user-defined check for line search step validity */
828: if (neP->postcheckstep) {
829: (*neP->postcheckstep)(snes,x,y,w,neP->postcheck,&changed_y,&changed_w);
830: if (changed_y) {
831: VecWAXPY(w,-1.0,y,x);
832: }
833: if (changed_y || changed_w) { /* recompute the function if the step has changed */
834: SNESComputeFunction(snes,w,g);
835: if (PetscExceptionValue(ierr)) {
837: }
838:
839: VecNormBegin(g,NORM_2,gnorm);
840: VecNormBegin(w,NORM_2,ynorm);
841: VecNormEnd(g,NORM_2,gnorm);
842: VecNormEnd(w,NORM_2,ynorm);
843: if (*gnorm != *gnorm) SETERRQ(PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
844: }
845: }
847: return(0);
848: }
850: /* -------------------------------------------------------------------------- */
853: /*@C
854: SNESLineSearchSet - Sets the line search routine to be used
855: by the method SNESLS.
857: Input Parameters:
858: + snes - nonlinear context obtained from SNESCreate()
859: . lsctx - optional user-defined context for use by line search
860: - func - pointer to int function
862: Collective on SNES
864: Available Routines:
865: + SNESLineSearchCubic() - default line search
866: . SNESLineSearchQuadratic() - quadratic line search
867: . SNESLineSearchNo() - the full Newton step (actually not a line search)
868: - SNESLineSearchNoNorms() - the full Newton step (calculating no norms; faster in parallel)
870: Options Database Keys:
871: + -snes_ls [cubic,quadratic,basic,basicnonorms] - Selects line search
872: . -snes_ls_alpha <alpha> - Sets alpha
873: . -snes_ls_maxstep <max> - Sets maxstep
874: - -snes_ls_steptol <steptol> - Sets steptol, this is the minimum step size that the line search code
875: will accept; min p[i]/x[i] < steptol. The -snes_stol <stol> is the minimum step length
876: the default convergence test will use and is based on 2-norm(p) < stol*2-norm(x)
878: Calling sequence of func:
879: .vb
880: func (SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,
881: PetscReal fnorm,PetscReal *ynorm,PetscReal *gnorm,PetscTruth *flag)
882: .ve
884: Input parameters for func:
885: + snes - nonlinear context
886: . lsctx - optional user-defined context for line search
887: . x - current iterate
888: . f - residual evaluated at x
889: . y - search direction
890: . w - work vector
891: - fnorm - 2-norm of f
893: Output parameters for func:
894: + g - residual evaluated at new iterate y
895: . w - new iterate
896: . gnorm - 2-norm of g
897: . ynorm - 2-norm of search length
898: - flag - set to PETSC_TRUE if the line search succeeds; PETSC_FALSE on failure.
900: Level: advanced
902: .keywords: SNES, nonlinear, set, line search, routine
904: .seealso: SNESLineSearchCubic(), SNESLineSearchQuadratic(), SNESLineSearchNo(), SNESLineSearchNoNorms(),
905: SNESLineSearchSetPostCheck(), SNESLineSearchSetParams(), SNESLineSearchGetParams(), SNESLineSearchSetPreCheck()
906: @*/
907: PetscErrorCode SNESLineSearchSet(SNES snes,PetscErrorCode (*func)(SNES,void*,Vec,Vec,Vec,Vec,Vec,PetscReal,PetscReal*,PetscReal*,PetscTruth*),void *lsctx)
908: {
909: PetscErrorCode ierr,(*f)(SNES,PetscErrorCode (*)(SNES,void*,Vec,Vec,Vec,Vec,Vec,PetscReal,PetscReal*,PetscReal*,PetscTruth*),void*);
912: PetscObjectQueryFunction((PetscObject)snes,"SNESLineSearchSet_C",(void (**)(void))&f);
913: if (f) {
914: (*f)(snes,func,lsctx);
915: }
916: return(0);
917: }
920: /* -------------------------------------------------------------------------- */
924: PetscErrorCode SNESLineSearchSet_LS(SNES snes,FCN2 func,void *lsctx)
925: {
927: ((SNES_LS *)(snes->data))->LineSearch = func;
928: ((SNES_LS *)(snes->data))->lsP = lsctx;
929: return(0);
930: }
932: /* -------------------------------------------------------------------------- */
935: /*@C
936: SNESLineSearchSetPostCheck - Sets a routine to check the validity of new iterate computed
937: by the line search routine in the Newton-based method SNESLS.
939: Input Parameters:
940: + snes - nonlinear context obtained from SNESCreate()
941: . func - pointer to function
942: - checkctx - optional user-defined context for use by step checking routine
944: Collective on SNES
946: Calling sequence of func:
947: .vb
948: int func (SNES snes, Vec x,Vec y,Vec w,void *checkctx, PetscTruth *changed_y,PetscTruth *changed_w)
949: .ve
950: where func returns an error code of 0 on success and a nonzero
951: on failure.
953: Input parameters for func:
954: + snes - nonlinear context
955: . checkctx - optional user-defined context for use by step checking routine
956: . x - previous iterate
957: . y - new search direction and length
958: - w - current candidate iterate
960: Output parameters for func:
961: + y - search direction (possibly changed)
962: . w - current iterate (possibly modified)
963: . changed_y - indicates search direction was changed by this routine
964: - changed_w - indicates current iterate was changed by this routine
966: Level: advanced
968: Notes: All line searches accept the new iterate computed by the line search checking routine.
970: Only one of changed_y and changed_w can be PETSC_TRUE
972: On input w = x + y
974: SNESLineSearchNo() and SNESLineSearchNoNorms() (1) compute a candidate iterate u_{i+1}, (2) pass control
975: to the checking routine, and then (3) compute the corresponding nonlinear
976: function f(u_{i+1}) with the (possibly altered) iterate u_{i+1}.
978: SNESLineSearchQuadratic() and SNESLineSearchCubic() (1) compute a candidate iterate u_{i+1} as well as a
979: candidate nonlinear function f(u_{i+1}), (2) pass control to the checking
980: routine, and then (3) force a re-evaluation of f(u_{i+1}) if any changes
981: were made to the candidate iterate in the checking routine (as indicated
982: by flag=PETSC_TRUE). The overhead of this extra function re-evaluation can be
983: very costly, so use this feature with caution!
985: .keywords: SNES, nonlinear, set, line search check, step check, routine
987: .seealso: SNESLineSearchSet(), SNESLineSearchSetPreCheck()
988: @*/
989: PetscErrorCode SNESLineSearchSetPostCheck(SNES snes,PetscErrorCode (*func)(SNES,Vec,Vec,Vec,void*,PetscTruth*,PetscTruth*),void *checkctx)
990: {
991: PetscErrorCode ierr,(*f)(SNES,PetscErrorCode (*)(SNES,Vec,Vec,Vec,void*,PetscTruth*,PetscTruth*),void*);
994: PetscObjectQueryFunction((PetscObject)snes,"SNESLineSearchSetPostCheck_C",(void (**)(void))&f);
995: if (f) {
996: (*f)(snes,func,checkctx);
997: }
998: return(0);
999: }
1002: /*@C
1003: SNESLineSearchSetPreCheck - Sets a routine to check the validity of a new direction given by the linear solve
1004: before the line search is called.
1006: Input Parameters:
1007: + snes - nonlinear context obtained from SNESCreate()
1008: . func - pointer to function
1009: - checkctx - optional user-defined context for use by step checking routine
1011: Collective on SNES
1013: Calling sequence of func:
1014: .vb
1015: int func (SNES snes, Vec x,Vec y,void *checkctx, PetscTruth *changed_y)
1016: .ve
1017: where func returns an error code of 0 on success and a nonzero
1018: on failure.
1020: Input parameters for func:
1021: + snes - nonlinear context
1022: . checkctx - optional user-defined context for use by step checking routine
1023: . x - previous iterate
1024: - y - new search direction and length
1026: Output parameters for func:
1027: + y - search direction (possibly changed)
1028: - changed_y - indicates search direction was changed by this routine
1030: Level: advanced
1032: Notes: All line searches accept the new iterate computed by the line search checking routine.
1034: .keywords: SNES, nonlinear, set, line search check, step check, routine
1036: .seealso: SNESLineSearchSet(), SNESLineSearchSetPostCheck(), SNESSetUpdate()
1037: @*/
1038: PetscErrorCode SNESLineSearchSetPreCheck(SNES snes,PetscErrorCode (*func)(SNES,Vec,Vec,void*,PetscTruth*),void *checkctx)
1039: {
1040: PetscErrorCode ierr,(*f)(SNES,PetscErrorCode (*)(SNES,Vec,Vec,void*,PetscTruth*),void*);
1043: PetscObjectQueryFunction((PetscObject)snes,"SNESLineSearchSetPreCheck_C",(void (**)(void))&f);
1044: if (f) {
1045: (*f)(snes,func,checkctx);
1046: }
1047: return(0);
1048: }
1050: /* -------------------------------------------------------------------------- */
1056: PetscErrorCode SNESLineSearchSetPostCheck_LS(SNES snes,FCN1 func,void *checkctx)
1057: {
1059: ((SNES_LS *)(snes->data))->postcheckstep = func;
1060: ((SNES_LS *)(snes->data))->postcheck = checkctx;
1061: return(0);
1062: }
1068: PetscErrorCode SNESLineSearchSetPreCheck_LS(SNES snes,FCN3 func,void *checkctx)
1069: {
1071: ((SNES_LS *)(snes->data))->precheckstep = func;
1072: ((SNES_LS *)(snes->data))->precheck = checkctx;
1073: return(0);
1074: }
1077: /*
1078: SNESView_LS - Prints info from the SNESLS data structure.
1080: Input Parameters:
1081: . SNES - the SNES context
1082: . viewer - visualization context
1084: Application Interface Routine: SNESView()
1085: */
1088: static PetscErrorCode SNESView_LS(SNES snes,PetscViewer viewer)
1089: {
1090: SNES_LS *ls = (SNES_LS *)snes->data;
1091: const char *cstr;
1093: PetscTruth iascii;
1096: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
1097: if (iascii) {
1098: if (ls->LineSearch == SNESLineSearchNo) cstr = "SNESLineSearchNo";
1099: else if (ls->LineSearch == SNESLineSearchQuadratic) cstr = "SNESLineSearchQuadratic";
1100: else if (ls->LineSearch == SNESLineSearchCubic) cstr = "SNESLineSearchCubic";
1101: else cstr = "unknown";
1102: PetscViewerASCIIPrintf(viewer," line search variant: %s\n",cstr);
1103: PetscViewerASCIIPrintf(viewer," alpha=%G, maxstep=%G, steptol=%G\n",ls->alpha,ls->maxstep,ls->steptol);
1104: } else {
1105: SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for SNES EQ LS",((PetscObject)viewer)->type_name);
1106: }
1107: return(0);
1108: }
1109: /* -------------------------------------------------------------------------- */
1110: /*
1111: SNESSetFromOptions_LS - Sets various parameters for the SNESLS method.
1113: Input Parameter:
1114: . snes - the SNES context
1116: Application Interface Routine: SNESSetFromOptions()
1117: */
1120: static PetscErrorCode SNESSetFromOptions_LS(SNES snes)
1121: {
1122: SNES_LS *ls = (SNES_LS *)snes->data;
1123: const char *lses[] = {"basic","basicnonorms","quadratic","cubic"};
1125: PetscInt indx;
1126: PetscTruth flg;
1129: PetscOptionsHead("SNES Line search options");
1130: PetscOptionsReal("-snes_ls_alpha","Function norm must decrease by","None",ls->alpha,&ls->alpha,0);
1131: PetscOptionsReal("-snes_ls_maxstep","Step must be less than","None",ls->maxstep,&ls->maxstep,0);
1132: PetscOptionsReal("-snes_ls_steptol","Step must be greater than","None",ls->steptol,&ls->steptol,0);
1134: PetscOptionsEList("-snes_ls","Line search used","SNESLineSearchSet",lses,4,"cubic",&indx,&flg);
1135: if (flg) {
1136: switch (indx) {
1137: case 0:
1138: SNESLineSearchSet(snes,SNESLineSearchNo,PETSC_NULL);
1139: break;
1140: case 1:
1141: SNESLineSearchSet(snes,SNESLineSearchNoNorms,PETSC_NULL);
1142: break;
1143: case 2:
1144: SNESLineSearchSet(snes,SNESLineSearchQuadratic,PETSC_NULL);
1145: break;
1146: case 3:
1147: SNESLineSearchSet(snes,SNESLineSearchCubic,PETSC_NULL);
1148: break;
1149: }
1150: }
1151: PetscOptionsTail();
1152: return(0);
1153: }
1154: /* -------------------------------------------------------------------------- */
1155: /*MC
1156: SNESLS - Newton based nonlinear solver that uses a line search
1158: Options Database:
1159: + -snes_ls [cubic,quadratic,basic,basicnonorms] - Selects line search
1160: . -snes_ls_alpha <alpha> - Sets alpha
1161: . -snes_ls_maxstep <max> - Sets maxstep
1162: - -snes_ls_steptol <steptol> - Sets steptol, this is the minimum step size that the line search code
1163: will accept; min p[i]/x[i] < steptol. The -snes_stol <stol> is the minimum step length
1164: the default convergence test will use and is based on 2-norm(p) < stol*2-norm(x)
1166: Notes: This is the default nonlinear solver in SNES
1168: Level: beginner
1170: .seealso: SNESCreate(), SNES, SNESSetType(), SNESTR, SNESLineSearchSet(),
1171: SNESLineSearchSetPostCheck(), SNESLineSearchNo(), SNESLineSearchCubic(), SNESLineSearchQuadratic(),
1172: SNESLineSearchSet(), SNESLineSearchNoNorms(), SNESLineSearchSetPreCheck()
1174: M*/
1178: PetscErrorCode SNESCreate_LS(SNES snes)
1179: {
1181: SNES_LS *neP;
1184: snes->ops->setup = SNESSetUp_LS;
1185: snes->ops->solve = SNESSolve_LS;
1186: snes->ops->destroy = SNESDestroy_LS;
1187: snes->ops->converged = SNESConverged_LS;
1188: snes->ops->setfromoptions = SNESSetFromOptions_LS;
1189: snes->ops->view = SNESView_LS;
1190: snes->nwork = 0;
1192: PetscNew(SNES_LS,&neP);
1193: PetscLogObjectMemory(snes,sizeof(SNES_LS));
1194: snes->data = (void*)neP;
1195: neP->alpha = 1.e-4;
1196: neP->maxstep = 1.e8;
1197: neP->steptol = 1.e-12;
1198: neP->LineSearch = SNESLineSearchCubic;
1199: neP->lsP = PETSC_NULL;
1200: neP->postcheckstep = PETSC_NULL;
1201: neP->postcheck = PETSC_NULL;
1202: neP->precheckstep = PETSC_NULL;
1203: neP->precheck = PETSC_NULL;
1205: PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSet_C","SNESLineSearchSet_LS",SNESLineSearchSet_LS);
1206: PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetPostCheck_C","SNESLineSearchSetPostCheck_LS",SNESLineSearchSetPostCheck_LS);
1207: PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetPreCheck_C","SNESLineSearchSetPreCheck_LS",SNESLineSearchSetPreCheck_LS);
1209: return(0);
1210: }