Actual source code: ls.c
1: /*$Id: ls.c,v 1.172 2001/08/07 03:04:11 balay Exp $*/
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: */
11: #undef __FUNCT__
13: int SNESLSCheckLocalMin_Private(Mat A,Vec F,Vec W,PetscReal fnorm,PetscTruth *ismin)
14: {
15: PetscReal a1;
16: int ierr;
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: PetscLogInfo(0,"SNESSolve_LS: || J^T F|| %g near zero implies found a local minimumn",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(PETSC_NULL,W);
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: PetscLogInfo(0,"SNESSolve_LS: (F^T J random)/(|| F ||*||random|| %g near zero implies found a local minimumn",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: */
49: #undef __FUNCT__
51: int SNESLSCheckResidual_Private(Mat A,Vec F,Vec X,Vec W1,Vec W2)
52: {
53: PetscReal a1,a2;
54: int ierr;
55: PetscTruth hastranspose;
56: PetscScalar mone = -1.0;
59: MatHasOperation(A,MATOP_MULT_TRANSPOSE,&hastranspose);
60: if (hastranspose) {
61: MatMult(A,X,W1);
62: VecAXPY(&mone,F,W1);
64: /* Compute || J^T W|| */
65: MatMultTranspose(A,W1,W2);
66: VecNorm(W1,NORM_2,&a1);
67: VecNorm(W2,NORM_2,&a2);
68: if (a1 != 0) {
69: PetscLogInfo(0,"SNESSolve_LS: ||J^T(F-Ax)||/||F-AX|| %g near zero implies inconsistent rhsn",a2/a1);
70: }
71: }
72: return(0);
73: }
75: /* --------------------------------------------------------------------
77: This file implements a truncated Newton method with a line search,
78: for solving a system of nonlinear equations, using the SLES, Vec,
79: and Mat interfaces for linear solvers, vectors, and matrices,
80: respectively.
82: The following basic routines are required for each nonlinear solver:
83: SNESCreate_XXX() - Creates a nonlinear solver context
84: SNESSetFromOptions_XXX() - Sets runtime options
85: SNESSolve_XXX() - Solves the nonlinear system
86: SNESDestroy_XXX() - Destroys the nonlinear solver context
87: The suffix "_XXX" denotes a particular implementation, in this case
88: we use _LS (e.g., SNESCreate_LS, SNESSolve_LS) for solving
89: systems of nonlinear equations with a line search (LS) method.
90: These routines are actually called via the common user interface
91: routines SNESCreate(), SNESSetFromOptions(), SNESSolve(), and
92: SNESDestroy(), so the application code interface remains identical
93: for all nonlinear solvers.
95: Another key routine is:
96: SNESSetUp_XXX() - Prepares for the use of a nonlinear solver
97: by setting data structures and options. The interface routine SNESSetUp()
98: is not usually called directly by the user, but instead is called by
99: SNESSolve() if necessary.
101: Additional basic routines are:
102: SNESView_XXX() - Prints details of runtime options that
103: have actually been used.
104: These are called by application codes via the interface routines
105: SNESView().
107: The various types of solvers (preconditioners, Krylov subspace methods,
108: nonlinear solvers, timesteppers) are all organized similarly, so the
109: above description applies to these categories also.
111: -------------------------------------------------------------------- */
112: /*
113: SNESSolve_LS - Solves a nonlinear system with a truncated Newton
114: method with a line search.
116: Input Parameters:
117: . snes - the SNES context
119: Output Parameter:
120: . outits - number of iterations until termination
122: Application Interface Routine: SNESSolve()
124: Notes:
125: This implements essentially a truncated Newton method with a
126: line search. By default a cubic backtracking line search
127: is employed, as described in the text "Numerical Methods for
128: Unconstrained Optimization and Nonlinear Equations" by Dennis
129: and Schnabel.
130: */
131: #undef __FUNCT__
133: int SNESSolve_LS(SNES snes,int *outits)
134: {
135: SNES_LS *neP = (SNES_LS*)snes->data;
136: int maxits,i,ierr,lits,lsfail;
137: MatStructure flg = DIFFERENT_NONZERO_PATTERN;
138: PetscReal fnorm,gnorm,xnorm,ynorm;
139: Vec Y,X,F,G,W,TMP;
142: snes->reason = SNES_CONVERGED_ITERATING;
144: maxits = snes->max_its; /* maximum number of iterations */
145: X = snes->vec_sol; /* solution vector */
146: F = snes->vec_func; /* residual vector */
147: Y = snes->work[0]; /* work vectors */
148: G = snes->work[1];
149: W = snes->work[2];
151: PetscObjectTakeAccess(snes);
152: snes->iter = 0;
153: PetscObjectGrantAccess(snes);
154: SNESComputeFunction(snes,X,F); /* F(X) */
155: VecNorm(F,NORM_2,&fnorm); /* fnorm <- ||F|| */
156: PetscObjectTakeAccess(snes);
157: snes->norm = fnorm;
158: PetscObjectGrantAccess(snes);
159: SNESLogConvHistory(snes,fnorm,0);
160: SNESMonitor(snes,0,fnorm);
162: if (fnorm < snes->atol) {*outits = 0; snes->reason = SNES_CONVERGED_FNORM_ABS; return(0);}
164: /* set parameter for default relative tolerance convergence test */
165: snes->ttol = fnorm*snes->rtol;
167: for (i=0; i<maxits; i++) {
169: /* Call general purpose update function */
170: if (snes->update != PETSC_NULL) {
171: (*snes->update)(snes, snes->iter);
172: }
174: /* Solve J Y = F, where J is Jacobian matrix */
175: SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);
176: SLESSetOperators(snes->sles,snes->jacobian,snes->jacobian_pre,flg);
177: SLESSolve(snes->sles,F,Y,&lits);
179: if (PetscLogPrintInfo){
180: SNESLSCheckResidual_Private(snes->jacobian,F,Y,G,W);
181: }
183: /* should check what happened to the linear solve? */
184: snes->linear_its += lits;
185: PetscLogInfo(snes,"SNESSolve_LS: iter=%d, linear solve iterations=%dn",snes->iter,lits);
187: /* Compute a (scaled) negative update in the line search routine:
188: Y <- X - lambda*Y
189: and evaluate G(Y) = function(Y))
190: */
191: VecCopy(Y,snes->vec_sol_update_always);
192: (*neP->LineSearch)(snes,neP->lsP,X,F,G,Y,W,fnorm,&ynorm,&gnorm,&lsfail);
193: PetscLogInfo(snes,"SNESSolve_LS: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lsfail=%dn",fnorm,gnorm,ynorm,lsfail);
195: TMP = F; F = G; snes->vec_func_always = F; G = TMP;
196: TMP = X; X = Y; snes->vec_sol_always = X; Y = TMP;
197: fnorm = gnorm;
199: PetscObjectTakeAccess(snes);
200: snes->iter = i+1;
201: snes->norm = fnorm;
202: PetscObjectGrantAccess(snes);
203: SNESLogConvHistory(snes,fnorm,lits);
204: SNESMonitor(snes,i+1,fnorm);
206: if (lsfail) {
207: PetscTruth ismin;
209: if (++snes->numFailures >= snes->maxFailures) {
210: snes->reason = SNES_DIVERGED_LS_FAILURE;
211: SNESLSCheckLocalMin_Private(snes->jacobian,F,W,fnorm,&ismin);
212: if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
213: break;
214: }
215: }
217: /* Test for convergence */
218: if (snes->converged) {
219: VecNorm(X,NORM_2,&xnorm); /* xnorm = || X || */
220: (*snes->converged)(snes,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);
221: if (snes->reason) {
222: break;
223: }
224: }
225: }
226: if (X != snes->vec_sol) {
227: VecCopy(X,snes->vec_sol);
228: }
229: if (F != snes->vec_func) {
230: VecCopy(F,snes->vec_func);
231: }
232: snes->vec_sol_always = snes->vec_sol;
233: snes->vec_func_always = snes->vec_func;
234: if (i == maxits) {
235: PetscLogInfo(snes,"SNESSolve_LS: Maximum number of iterations has been reached: %dn",maxits);
236: i--;
237: snes->reason = SNES_DIVERGED_MAX_IT;
238: }
239: PetscObjectTakeAccess(snes);
240: PetscObjectGrantAccess(snes);
241: *outits = i+1;
242: return(0);
243: }
244: /* -------------------------------------------------------------------------- */
245: /*
246: SNESSetUp_LS - Sets up the internal data structures for the later use
247: of the SNESLS nonlinear solver.
249: Input Parameter:
250: . snes - the SNES context
251: . x - the solution vector
253: Application Interface Routine: SNESSetUp()
255: Notes:
256: For basic use of the SNES solvers, the user need not explicitly call
257: SNESSetUp(), since these actions will automatically occur during
258: the call to SNESSolve().
259: */
260: #undef __FUNCT__
262: int SNESSetUp_LS(SNES snes)
263: {
267: snes->nwork = 4;
268: VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work);
269: PetscLogObjectParents(snes,snes->nwork,snes->work);
270: snes->vec_sol_update_always = snes->work[3];
271: return(0);
272: }
273: /* -------------------------------------------------------------------------- */
274: /*
275: SNESDestroy_LS - Destroys the private SNES_LS context that was created
276: with SNESCreate_LS().
278: Input Parameter:
279: . snes - the SNES context
281: Application Interface Routine: SNESDestroy()
282: */
283: #undef __FUNCT__
285: int SNESDestroy_LS(SNES snes)
286: {
287: int ierr;
290: if (snes->nwork) {
291: VecDestroyVecs(snes->work,snes->nwork);
292: }
293: PetscFree(snes->data);
294: return(0);
295: }
296: /* -------------------------------------------------------------------------- */
297: #undef __FUNCT__
300: /*@C
301: SNESNoLineSearch - This routine is not a line search at all;
302: it simply uses the full Newton step. Thus, this routine is intended
303: to serve as a template and is not recommended for general use.
305: Collective on SNES and Vec
307: Input Parameters:
308: + snes - nonlinear context
309: . lsctx - optional context for line search (not used here)
310: . x - current iterate
311: . f - residual evaluated at x
312: . y - search direction (contains new iterate on output)
313: . w - work vector
314: - fnorm - 2-norm of f
316: Output Parameters:
317: + g - residual evaluated at new iterate y
318: . y - new iterate (contains search direction on input)
319: . gnorm - 2-norm of g
320: . ynorm - 2-norm of search length
321: - flag - set to 0, indicating a successful line search
323: Options Database Key:
324: . -snes_ls basic - Activates SNESNoLineSearch()
326: Level: advanced
328: .keywords: SNES, nonlinear, line search, cubic
330: .seealso: SNESCubicLineSearch(), SNESQuadraticLineSearch(),
331: SNESSetLineSearch(), SNESNoLineSearchNoNorms()
332: @*/
333: int SNESNoLineSearch(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal *ynorm,PetscReal *gnorm,int *flag)
334: {
335: int ierr;
336: PetscScalar mone = -1.0;
337: SNES_LS *neP = (SNES_LS*)snes->data;
338: PetscTruth change_y = PETSC_FALSE;
341: *flag = 0;
342: PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);
343: VecNorm(y,NORM_2,ynorm); /* ynorm = || y || */
344: VecAYPX(&mone,x,y); /* y <- y - x */
345: if (neP->CheckStep) {
346: (*neP->CheckStep)(snes,neP->checkP,y,&change_y);
347: }
348: SNESComputeFunction(snes,y,g); /* Compute F(y) */
349: VecNorm(g,NORM_2,gnorm); /* gnorm = || g || */
350: PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);
351: return(0);
352: }
353: /* -------------------------------------------------------------------------- */
354: #undef __FUNCT__
357: /*@C
358: SNESNoLineSearchNoNorms - This routine is not a line search at
359: all; it simply uses the full Newton step. This version does not
360: even compute the norm of the function or search direction; this
361: is intended only when you know the full step is fine and are
362: not checking for convergence of the nonlinear iteration (for
363: example, you are running always for a fixed number of Newton steps).
365: Collective on SNES and Vec
367: Input Parameters:
368: + snes - nonlinear context
369: . lsctx - optional context for line search (not used here)
370: . x - current iterate
371: . f - residual evaluated at x
372: . y - search direction (contains new iterate on output)
373: . w - work vector
374: - fnorm - 2-norm of f
376: Output Parameters:
377: + g - residual evaluated at new iterate y
378: . gnorm - not changed
379: . ynorm - not changed
380: - flag - set to 0, indicating a successful line search
382: Options Database Key:
383: . -snes_ls basicnonorms - Activates SNESNoLineSearchNoNorms()
385: Notes:
386: SNESNoLineSearchNoNorms() must be used in conjunction with
387: either the options
388: $ -snes_no_convergence_test -snes_max_it <its>
389: or alternatively a user-defined custom test set via
390: SNESSetConvergenceTest(); or a -snes_max_it of 1,
391: otherwise, the SNES solver will generate an error.
393: During the final iteration this will not evaluate the function at
394: the solution point. This is to save a function evaluation while
395: using pseudo-timestepping.
397: The residual norms printed by monitoring routines such as
398: SNESDefaultMonitor() (as activated via -snes_monitor) will not be
399: correct, since they are not computed.
401: Level: advanced
403: .keywords: SNES, nonlinear, line search, cubic
405: .seealso: SNESCubicLineSearch(), SNESQuadraticLineSearch(),
406: SNESSetLineSearch(), SNESNoLineSearch()
407: @*/
408: int SNESNoLineSearchNoNorms(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal *ynorm,PetscReal *gnorm,int *flag)
409: {
410: int ierr;
411: PetscScalar mone = -1.0;
412: SNES_LS *neP = (SNES_LS*)snes->data;
413: PetscTruth change_y = PETSC_FALSE;
416: *flag = 0;
417: PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);
418: VecAYPX(&mone,x,y); /* y <- y - x */
419: if (neP->CheckStep) {
420: (*neP->CheckStep)(snes,neP->checkP,y,&change_y);
421: }
422:
423: /* don't evaluate function the last time through */
424: if (snes->iter < snes->max_its-1) {
425: SNESComputeFunction(snes,y,g); /* Compute F(y) */
426: }
427: PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);
428: return(0);
429: }
430: /* -------------------------------------------------------------------------- */
431: #undef __FUNCT__
433: /*@C
434: SNESCubicLineSearch - Performs a cubic line search (default line search method).
436: Collective on SNES
438: Input Parameters:
439: + snes - nonlinear context
440: . lsctx - optional context for line search (not used here)
441: . x - current iterate
442: . f - residual evaluated at x
443: . y - search direction (contains new iterate on output)
444: . w - work vector
445: - fnorm - 2-norm of f
447: Output Parameters:
448: + g - residual evaluated at new iterate y
449: . y - new iterate (contains search direction on input)
450: . gnorm - 2-norm of g
451: . ynorm - 2-norm of search length
452: - flag - 0 if line search succeeds; -1 on failure.
454: Options Database Key:
455: $ -snes_ls cubic - Activates SNESCubicLineSearch()
457: Notes:
458: This line search is taken from "Numerical Methods for Unconstrained
459: Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325.
461: Level: advanced
463: .keywords: SNES, nonlinear, line search, cubic
465: .seealso: SNESQuadraticLineSearch(), SNESNoLineSearch(), SNESSetLineSearch(), SNESNoLineSearchNoNorms()
466: @*/
467: int SNESCubicLineSearch(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal *ynorm,PetscReal *gnorm,int *flag)
468: {
469: /*
470: Note that for line search purposes we work with with the related
471: minimization problem:
472: min z(x): R^n -> R,
473: where z(x) = .5 * fnorm*fnorm, and fnorm = || f ||_2.
474: */
475:
476: PetscReal steptol,initslope,lambdaprev,gnormprev,a,b,d,t1,t2,rellength;
477: PetscReal maxstep,minlambda,alpha,lambda,lambdatemp,lambdaneg;
478: #if defined(PETSC_USE_COMPLEX)
479: PetscScalar cinitslope,clambda;
480: #endif
481: int ierr,count;
482: SNES_LS *neP = (SNES_LS*)snes->data;
483: PetscScalar mone = -1.0,scale;
484: PetscTruth change_y = PETSC_FALSE;
487: PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);
488: *flag = 0;
489: alpha = neP->alpha;
490: maxstep = neP->maxstep;
491: steptol = neP->steptol;
493: VecNorm(y,NORM_2,ynorm);
494: if (*ynorm < snes->atol) {
495: PetscLogInfo(snes,"SNESCubicLineSearch: Search direction and size are nearly 0n");
496: *gnorm = fnorm;
497: VecCopy(x,y);
498: VecCopy(f,g);
499: goto theend1;
500: }
501: if (*ynorm > maxstep) { /* Step too big, so scale back */
502: scale = maxstep/(*ynorm);
503: #if defined(PETSC_USE_COMPLEX)
504: PetscLogInfo(snes,"SNESCubicLineSearch: Scaling step by %g old ynorm %gn",PetscRealPart(scale),*ynorm);
505: #else
506: PetscLogInfo(snes,"SNESCubicLineSearch: Scaling step by %g old ynorm %gn",scale,*ynorm);
507: #endif
508: VecScale(&scale,y);
509: *ynorm = maxstep;
510: }
511: ierr = VecMaxPointwiseDivide(y,x,&rellength);
512: minlambda = steptol/rellength;
513: MatMult(snes->jacobian,y,w);
514: #if defined(PETSC_USE_COMPLEX)
515: VecDot(f,w,&cinitslope);
516: initslope = PetscRealPart(cinitslope);
517: #else
518: VecDot(f,w,&initslope);
519: #endif
520: if (initslope > 0.0) initslope = -initslope;
521: if (initslope == 0.0) initslope = -1.0;
523: VecCopy(y,w);
524: VecAYPX(&mone,x,w);
525: SNESComputeFunction(snes,w,g);
526: VecNorm(g,NORM_2,gnorm);
527: if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + alpha*initslope) { /* Sufficient reduction */
528: VecCopy(w,y);
529: PetscLogInfo(snes,"SNESCubicLineSearch: Using full stepn");
530: goto theend1;
531: }
533: /* Fit points with quadratic */
534: lambda = 1.0;
535: lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope);
536: lambdaprev = lambda;
537: gnormprev = *gnorm;
538: if (lambdatemp > .5*lambda) lambdatemp = .5*lambda;
539: if (lambdatemp <= .1*lambda) lambda = .1*lambda;
540: else lambda = lambdatemp;
541: ierr = VecCopy(x,w);
542: lambdaneg = -lambda;
543: #if defined(PETSC_USE_COMPLEX)
544: clambda = lambdaneg; VecAXPY(&clambda,y,w);
545: #else
546: VecAXPY(&lambdaneg,y,w);
547: #endif
548: SNESComputeFunction(snes,w,g);
549: VecNorm(g,NORM_2,gnorm);
550: if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*alpha*initslope) { /* sufficient reduction */
551: VecCopy(w,y);
552: PetscLogInfo(snes,"SNESCubicLineSearch: Quadratically determined step, lambda=%18.16en",lambda);
553: goto theend1;
554: }
556: /* Fit points with cubic */
557: count = 1;
558: while (PETSC_TRUE) {
559: if (lambda <= minlambda) { /* bad luck; use full step */
560: PetscLogInfo(snes,"SNESCubicLineSearch:Unable to find good step length! %d n",count);
561: PetscLogInfo(snes,"SNESCubicLineSearch:fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16en",fnorm,*gnorm,*ynorm,lambda,initslope);
562: VecCopy(x,y);
563: *flag = -1; break;
564: }
565: t1 = .5*((*gnorm)*(*gnorm) - fnorm*fnorm) - lambda*initslope;
566: t2 = .5*(gnormprev*gnormprev - fnorm*fnorm) - lambdaprev*initslope;
567: a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
568: b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
569: d = b*b - 3*a*initslope;
570: if (d < 0.0) d = 0.0;
571: if (a == 0.0) {
572: lambdatemp = -initslope/(2.0*b);
573: } else {
574: lambdatemp = (-b + sqrt(d))/(3.0*a);
575: }
576: lambdaprev = lambda;
577: gnormprev = *gnorm;
578: if (lambdatemp > .5*lambda) lambdatemp = .5*lambda;
579: if (lambdatemp <= .1*lambda) lambda = .1*lambda;
580: else lambda = lambdatemp;
581: VecCopy(x,w);
582: lambdaneg = -lambda;
583: #if defined(PETSC_USE_COMPLEX)
584: clambda = lambdaneg;
585: VecAXPY(&clambda,y,w);
586: #else
587: VecAXPY(&lambdaneg,y,w);
588: #endif
589: SNESComputeFunction(snes,w,g);
590: VecNorm(g,NORM_2,gnorm);
591: if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*alpha*initslope) { /* is reduction enough? */
592: VecCopy(w,y);
593: PetscLogInfo(snes,"SNESCubicLineSearch: Cubically determined step, lambda=%18.16en",lambda);
594: goto theend1;
595: } else {
596: PetscLogInfo(snes,"SNESCubicLineSearch: Cubic step no good, shrinking lambda, lambda=%18.16en",lambda);
597: }
598: count++;
599: }
600: theend1:
601: /* Optional user-defined check for line search step validity */
602: if (neP->CheckStep) {
603: (*neP->CheckStep)(snes,neP->checkP,y,&change_y);
604: if (change_y == PETSC_TRUE) { /* recompute the function if the step has changed */
605: SNESComputeFunction(snes,y,g);
606: VecNormBegin(y,NORM_2,ynorm);
607: VecNormBegin(g,NORM_2,gnorm);
608: VecNormEnd(y,NORM_2,ynorm);
609: VecNormEnd(g,NORM_2,gnorm);
610: }
611: }
612: PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);
613: return(0);
614: }
615: /* -------------------------------------------------------------------------- */
616: #undef __FUNCT__
618: /*@C
619: SNESQuadraticLineSearch - Performs a quadratic line search.
621: Collective on SNES and Vec
623: Input Parameters:
624: + snes - the SNES context
625: . lsctx - optional context for line search (not used here)
626: . x - current iterate
627: . f - residual evaluated at x
628: . y - search direction (contains new iterate on output)
629: . w - work vector
630: - fnorm - 2-norm of f
632: Output Parameters:
633: + g - residual evaluated at new iterate y
634: . y - new iterate (contains search direction on input)
635: . gnorm - 2-norm of g
636: . ynorm - 2-norm of search length
637: - flag - 0 if line search succeeds; -1 on failure.
639: Options Database Key:
640: . -snes_ls quadratic - Activates SNESQuadraticLineSearch()
642: Notes:
643: Use SNESSetLineSearch() to set this routine within the SNESLS method.
645: Level: advanced
647: .keywords: SNES, nonlinear, quadratic, line search
649: .seealso: SNESCubicLineSearch(), SNESNoLineSearch(), SNESSetLineSearch(), SNESNoLineSearchNoNorms()
650: @*/
651: int SNESQuadraticLineSearch(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal *ynorm,PetscReal *gnorm,int *flag)
652: {
653: /*
654: Note that for line search purposes we work with with the related
655: minimization problem:
656: min z(x): R^n -> R,
657: where z(x) = .5 * fnorm*fnorm,and fnorm = || f ||_2.
658: */
659: PetscReal steptol,initslope,maxstep,minlambda,alpha,lambda,lambdatemp,lambdaneg,rellength;
660: #if defined(PETSC_USE_COMPLEX)
661: PetscScalar cinitslope,clambda;
662: #endif
663: int ierr,count;
664: SNES_LS *neP = (SNES_LS*)snes->data;
665: PetscScalar mone = -1.0,scale;
666: PetscTruth change_y = PETSC_FALSE;
669: PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);
670: *flag = 0;
671: alpha = neP->alpha;
672: maxstep = neP->maxstep;
673: steptol = neP->steptol;
675: VecNorm(y,NORM_2,ynorm);
676: if (*ynorm < snes->atol) {
677: PetscLogInfo(snes,"SNESQuadraticLineSearch: Search direction and size is 0n");
678: *gnorm = fnorm;
679: VecCopy(x,y);
680: VecCopy(f,g);
681: goto theend2;
682: }
683: if (*ynorm > maxstep) { /* Step too big, so scale back */
684: scale = maxstep/(*ynorm);
685: VecScale(&scale,y);
686: *ynorm = maxstep;
687: }
688: ierr = VecMaxPointwiseDivide(y,x,&rellength);
689: minlambda = steptol/rellength;
690: MatMult(snes->jacobian,y,w);
691: #if defined(PETSC_USE_COMPLEX)
692: VecDot(f,w,&cinitslope);
693: initslope = PetscRealPart(cinitslope);
694: #else
695: VecDot(f,w,&initslope);
696: #endif
697: if (initslope > 0.0) initslope = -initslope;
698: if (initslope == 0.0) initslope = -1.0;
700: VecCopy(y,w);
701: VecAYPX(&mone,x,w);
702: SNESComputeFunction(snes,w,g);
703: VecNorm(g,NORM_2,gnorm);
704: if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + alpha*initslope) { /* Sufficient reduction */
705: VecCopy(w,y);
706: PetscLogInfo(snes,"SNESQuadraticLineSearch: Using full stepn");
707: goto theend2;
708: }
710: /* Fit points with quadratic */
711: lambda = 1.0;
712: count = 1;
713: while (PETSC_TRUE) {
714: if (lambda <= minlambda) { /* bad luck; use full step */
715: PetscLogInfo(snes,"SNESQuadraticLineSearch:Unable to find good step length! %d n",count);
716: PetscLogInfo(snes,"SNESQuadraticLineSearch:fnorm=%g, gnorm=%g, ynorm=%g, lambda=%g, initial slope=%gn",fnorm,*gnorm,*ynorm,lambda,initslope);
717: VecCopy(x,y);
718: *flag = -1; break;
719: }
720: lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope);
721: if (lambdatemp > .5*lambda) lambdatemp = .5*lambda;
722: if (lambdatemp <= .1*lambda) lambda = .1*lambda;
723: else lambda = lambdatemp;
724: VecCopy(x,w);
725: lambdaneg = -lambda;
726: #if defined(PETSC_USE_COMPLEX)
727: clambda = lambdaneg; VecAXPY(&clambda,y,w);
728: #else
729: VecAXPY(&lambdaneg,y,w);
730: #endif
731: SNESComputeFunction(snes,w,g);
732: VecNorm(g,NORM_2,gnorm);
733: if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*alpha*initslope) { /* sufficient reduction */
734: VecCopy(w,y);
735: PetscLogInfo(snes,"SNESQuadraticLineSearch:Quadratically determined step, lambda=%gn",lambda);
736: break;
737: }
738: count++;
739: }
740: theend2:
741: /* Optional user-defined check for line search step validity */
742: if (neP->CheckStep) {
743: (*neP->CheckStep)(snes,neP->checkP,y,&change_y);
744: if (change_y == PETSC_TRUE) { /* recompute the function if the step has changed */
745: SNESComputeFunction(snes,y,g);
746: VecNormBegin(y,NORM_2,ynorm);
747: VecNormBegin(g,NORM_2,gnorm);
748: VecNormEnd(y,NORM_2,ynorm);
749: VecNormEnd(g,NORM_2,gnorm);
750: }
751: }
752: PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);
753: return(0);
754: }
755: /* -------------------------------------------------------------------------- */
756: #undef __FUNCT__
758: /*@C
759: SNESSetLineSearch - Sets the line search routine to be used
760: by the method SNESLS.
762: Input Parameters:
763: + snes - nonlinear context obtained from SNESCreate()
764: . lsctx - optional user-defined context for use by line search
765: - func - pointer to int function
767: Collective on SNES
769: Available Routines:
770: + SNESCubicLineSearch() - default line search
771: . SNESQuadraticLineSearch() - quadratic line search
772: . SNESNoLineSearch() - the full Newton step (actually not a line search)
773: - SNESNoLineSearchNoNorms() - the full Newton step (calculating no norms; faster in parallel)
775: Options Database Keys:
776: + -snes_ls [cubic,quadratic,basic,basicnonorms] - Selects line search
777: . -snes_ls_alpha <alpha> - Sets alpha
778: . -snes_ls_maxstep <max> - Sets maxstep
779: - -snes_ls_steptol <steptol> - Sets steptol, this is the minimum step size that the line search code
780: will accept; min p[i]/x[i] < steptol. The -snes_stol <stol> is the minimum step length
781: the default convergence test will use and is based on 2-norm(p) < stol*2-norm(x)
783: Calling sequence of func:
784: .vb
785: func (SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,
786: PetscReal fnorm,PetscReal *ynorm,PetscReal *gnorm,*flag)
787: .ve
789: Input parameters for func:
790: + snes - nonlinear context
791: . lsctx - optional user-defined context for line search
792: . x - current iterate
793: . f - residual evaluated at x
794: . y - search direction (contains new iterate on output)
795: . w - work vector
796: - fnorm - 2-norm of f
798: Output parameters for func:
799: + g - residual evaluated at new iterate y
800: . y - new iterate (contains search direction on input)
801: . gnorm - 2-norm of g
802: . ynorm - 2-norm of search length
803: - flag - set to 0 if the line search succeeds; a nonzero integer
804: on failure.
806: Level: advanced
808: .keywords: SNES, nonlinear, set, line search, routine
810: .seealso: SNESCubicLineSearch(), SNESQuadraticLineSearch(), SNESNoLineSearch(), SNESNoLineSearchNoNorms(),
811: SNESSetLineSearchCheck(), SNESSetLineSearchParams(), SNESGetLineSearchParams()
812: @*/
813: int SNESSetLineSearch(SNES snes,int (*func)(SNES,void*,Vec,Vec,Vec,Vec,Vec,PetscReal,PetscReal*,PetscReal*,int*),void *lsctx)
814: {
815: int ierr,(*f)(SNES,int (*)(SNES,void*,Vec,Vec,Vec,Vec,Vec,PetscReal,PetscReal*,PetscReal*,int*),void*);
818: PetscObjectQueryFunction((PetscObject)snes,"SNESSetLineSearch_C",(void (**)(void))&f);
819: if (f) {
820: (*f)(snes,func,lsctx);
821: }
822: return(0);
823: }
824: /* -------------------------------------------------------------------------- */
825: EXTERN_C_BEGIN
826: #undef __FUNCT__
828: int SNESSetLineSearch_LS(SNES snes,int (*func)(SNES,void*,Vec,Vec,Vec,Vec,Vec,
829: PetscReal,PetscReal*,PetscReal*,int*),void *lsctx)
830: {
832: ((SNES_LS *)(snes->data))->LineSearch = func;
833: ((SNES_LS *)(snes->data))->lsP = lsctx;
834: return(0);
835: }
836: EXTERN_C_END
837: /* -------------------------------------------------------------------------- */
838: #undef __FUNCT__
840: /*@C
841: SNESSetLineSearchCheck - Sets a routine to check the validity of new iterate computed
842: by the line search routine in the Newton-based method SNESLS.
844: Input Parameters:
845: + snes - nonlinear context obtained from SNESCreate()
846: . func - pointer to int function
847: - checkctx - optional user-defined context for use by step checking routine
849: Collective on SNES
851: Calling sequence of func:
852: .vb
853: int func (SNES snes, void *checkctx, Vec x, PetscTruth *flag)
854: .ve
855: where func returns an error code of 0 on success and a nonzero
856: on failure.
858: Input parameters for func:
859: + snes - nonlinear context
860: . checkctx - optional user-defined context for use by step checking routine
861: - x - current candidate iterate
863: Output parameters for func:
864: + x - current iterate (possibly modified)
865: - flag - flag indicating whether x has been modified (either
866: PETSC_TRUE of PETSC_FALSE)
868: Level: advanced
870: Notes:
871: SNESNoLineSearch() and SNESNoLineSearchNoNorms() accept the new
872: iterate computed by the line search checking routine. In particular,
873: these routines (1) compute a candidate iterate u_{i+1}, (2) pass control
874: to the checking routine, and then (3) compute the corresponding nonlinear
875: function f(u_{i+1}) with the (possibly altered) iterate u_{i+1}.
877: SNESQuadraticLineSearch() and SNESCubicLineSearch() also accept the
878: new iterate computed by the line search checking routine. In particular,
879: these routines (1) compute a candidate iterate u_{i+1} as well as a
880: candidate nonlinear function f(u_{i+1}), (2) pass control to the checking
881: routine, and then (3) force a re-evaluation of f(u_{i+1}) if any changes
882: were made to the candidate iterate in the checking routine (as indicated
883: by flag=PETSC_TRUE). The overhead of this function re-evaluation can be
884: very costly, so use this feature with caution!
886: .keywords: SNES, nonlinear, set, line search check, step check, routine
888: .seealso: SNESSetLineSearch()
889: @*/
890: int SNESSetLineSearchCheck(SNES snes,int (*func)(SNES,void*,Vec,PetscTruth*),void *checkctx)
891: {
892: int ierr,(*f)(SNES,int (*)(SNES,void*,Vec,PetscTruth*),void*);
895: PetscObjectQueryFunction((PetscObject)snes,"SNESSetLineSearchCheck_C",(void (**)(void))&f);
896: if (f) {
897: (*f)(snes,func,checkctx);
898: }
899: return(0);
900: }
901: /* -------------------------------------------------------------------------- */
902: EXTERN_C_BEGIN
903: #undef __FUNCT__
905: int SNESSetLineSearchCheck_LS(SNES snes,int (*func)(SNES,void*,Vec,PetscTruth*),void *checkctx)
906: {
908: ((SNES_LS *)(snes->data))->CheckStep = func;
909: ((SNES_LS *)(snes->data))->checkP = checkctx;
910: return(0);
911: }
912: EXTERN_C_END
913: /* -------------------------------------------------------------------------- */
914: /*
915: SNESPrintHelp_LS - Prints all options for the SNES_LS method.
917: Input Parameter:
918: . snes - the SNES context
920: Application Interface Routine: SNESPrintHelp()
921: */
922: #undef __FUNCT__
924: static int SNESPrintHelp_LS(SNES snes,char *p)
925: {
926: SNES_LS *ls = (SNES_LS *)snes->data;
929: (*PetscHelpPrintf)(snes->comm," method SNES_LS (ls) for systems of nonlinear equations:n");
930: (*PetscHelpPrintf)(snes->comm," %ssnes_ls [cubic,quadratic,basic,basicnonorms,...]n",p);
931: (*PetscHelpPrintf)(snes->comm," %ssnes_ls_alpha <alpha> (default %g)n",p,ls->alpha);
932: (*PetscHelpPrintf)(snes->comm," %ssnes_ls_maxstep <max> (default %g)n",p,ls->maxstep);
933: (*PetscHelpPrintf)(snes->comm," %ssnes_ls_steptol <tol> (default %g)n",p,ls->steptol);
934: return(0);
935: }
937: /*
938: SNESView_LS - Prints info from the SNESLS data structure.
940: Input Parameters:
941: . SNES - the SNES context
942: . viewer - visualization context
944: Application Interface Routine: SNESView()
945: */
946: #undef __FUNCT__
948: static int SNESView_LS(SNES snes,PetscViewer viewer)
949: {
950: SNES_LS *ls = (SNES_LS *)snes->data;
951: char *cstr;
952: int ierr;
953: PetscTruth isascii;
956: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
957: if (isascii) {
958: if (ls->LineSearch == SNESNoLineSearch) cstr = "SNESNoLineSearch";
959: else if (ls->LineSearch == SNESQuadraticLineSearch) cstr = "SNESQuadraticLineSearch";
960: else if (ls->LineSearch == SNESCubicLineSearch) cstr = "SNESCubicLineSearch";
961: else cstr = "unknown";
962: PetscViewerASCIIPrintf(viewer," line search variant: %sn",cstr);
963: PetscViewerASCIIPrintf(viewer," alpha=%g, maxstep=%g, steptol=%gn",ls->alpha,ls->maxstep,ls->steptol);
964: } else {
965: SETERRQ1(1,"Viewer type %s not supported for SNES EQ LS",((PetscObject)viewer)->type_name);
966: }
967: return(0);
968: }
969: /* -------------------------------------------------------------------------- */
970: /*
971: SNESSetFromOptions_LS - Sets various parameters for the SNESLS method.
973: Input Parameter:
974: . snes - the SNES context
976: Application Interface Routine: SNESSetFromOptions()
977: */
978: #undef __FUNCT__
980: static int SNESSetFromOptions_LS(SNES snes)
981: {
982: SNES_LS *ls = (SNES_LS *)snes->data;
983: char ver[16],*lses[] = {"basic","basicnonorms","quadratic","cubic"};
984: int ierr;
985: PetscTruth flg;
988: PetscOptionsHead("SNES Line search options");
989: PetscOptionsReal("-snes_ls_alpha","Function norm must decrease by","None",ls->alpha,&ls->alpha,0);
990: PetscOptionsReal("-snes_ls_maxstep","Step must be less than","None",ls->maxstep,&ls->maxstep,0);
991: PetscOptionsReal("-snes_ls_steptol","Step must be greater than","None",ls->steptol,&ls->steptol,0);
993: PetscOptionsEList("-snes_ls","Line search used","SNESSetLineSearch",lses,4,"cubic",ver,16,&flg);
994: if (flg) {
995: PetscTruth isbasic,isnonorms,isquad,iscubic;
997: PetscStrcmp(ver,lses[0],&isbasic);
998: PetscStrcmp(ver,lses[1],&isnonorms);
999: PetscStrcmp(ver,lses[2],&isquad);
1000: PetscStrcmp(ver,lses[3],&iscubic);
1002: if (isbasic) {
1003: SNESSetLineSearch(snes,SNESNoLineSearch,PETSC_NULL);
1004: } else if (isnonorms) {
1005: SNESSetLineSearch(snes,SNESNoLineSearchNoNorms,PETSC_NULL);
1006: } else if (isquad) {
1007: SNESSetLineSearch(snes,SNESQuadraticLineSearch,PETSC_NULL);
1008: } else if (iscubic) {
1009: SNESSetLineSearch(snes,SNESCubicLineSearch,PETSC_NULL);
1010: }
1011: else {SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Unknown line search");}
1012: }
1013: PetscOptionsTail();
1014: return(0);
1015: }
1016: /* -------------------------------------------------------------------------- */
1017: /*
1018: SNESCreate_LS - Creates a nonlinear solver context for the SNESLS method,
1019: SNES_LS, and sets this as the private data within the generic nonlinear solver
1020: context, SNES, that was created within SNESCreate().
1023: Input Parameter:
1024: . snes - the SNES context
1026: Application Interface Routine: SNESCreate()
1027: */
1028: EXTERN_C_BEGIN
1029: #undef __FUNCT__
1031: int SNESCreate_LS(SNES snes)
1032: {
1033: int ierr;
1034: SNES_LS *neP;
1037: snes->setup = SNESSetUp_LS;
1038: snes->solve = SNESSolve_LS;
1039: snes->destroy = SNESDestroy_LS;
1040: snes->converged = SNESConverged_LS;
1041: snes->printhelp = SNESPrintHelp_LS;
1042: snes->setfromoptions = SNESSetFromOptions_LS;
1043: snes->view = SNESView_LS;
1044: snes->nwork = 0;
1046: ierr = PetscNew(SNES_LS,&neP);
1047: PetscLogObjectMemory(snes,sizeof(SNES_LS));
1048: snes->data = (void*)neP;
1049: neP->alpha = 1.e-4;
1050: neP->maxstep = 1.e8;
1051: neP->steptol = 1.e-12;
1052: neP->LineSearch = SNESCubicLineSearch;
1053: neP->lsP = PETSC_NULL;
1054: neP->CheckStep = PETSC_NULL;
1055: neP->checkP = PETSC_NULL;
1057: PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESSetLineSearch_C","SNESSetLineSearch_LS",SNESSetLineSearch_LS);
1058: PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESSetLineSearchCheck_C","SNESSetLineSearchCheck_LS",SNESSetLineSearchCheck_LS);
1060: return(0);
1061: }
1062: EXTERN_C_END