Actual source code: linesearchcp.c
petsc-3.9.4 2018-09-11
1: #include <petsc/private/linesearchimpl.h>
2: #include <petscsnes.h>
4: static PetscErrorCode SNESLineSearchApply_CP(SNESLineSearch linesearch)
5: {
6: PetscBool changed_y, changed_w;
8: Vec X, Y, F, W;
9: SNES snes;
10: PetscReal xnorm, ynorm, gnorm, steptol, atol, rtol, ltol, maxstep;
12: PetscReal lambda, lambda_old, lambda_update, delLambda;
13: PetscScalar fty, fty_init, fty_old, fty_mid1, fty_mid2, s;
14: PetscInt i, max_its;
16: PetscViewer monitor;
19: SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, NULL);
20: SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);
21: SNESLineSearchGetSNES(linesearch, &snes);
22: SNESLineSearchGetLambda(linesearch, &lambda);
23: SNESLineSearchGetTolerances(linesearch, &steptol, &maxstep, &rtol, &atol, <ol, &max_its);
24: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED);
25: SNESLineSearchGetDefaultMonitor(linesearch, &monitor);
27: /* precheck */
28: SNESLineSearchPreCheck(linesearch,X,Y,&changed_y);
29: lambda_old = 0.0;
31: VecDot(F,Y,&fty_old);
32: if (PetscAbsScalar(fty_old) < atol) {
33: if (monitor) {
34: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
35: PetscViewerASCIIPrintf(monitor," Line search terminated ended at initial point because dot(F,Y) = %g < atol = %g\n",(double)PetscAbsScalar(fty_old), (double)atol);
36: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
37: }
38: return(0);
39: }
41: fty_init = fty_old;
43: for (i = 0; i < max_its; i++) {
44: /* compute the norm at lambda */
45: VecCopy(X, W);
46: VecAXPY(W, -lambda, Y);
47: if (linesearch->ops->viproject) {
48: (*linesearch->ops->viproject)(snes, W);
49: }
50: (*linesearch->ops->snesfunc)(snes,W,F);
51: VecDot(F,Y,&fty);
53: delLambda = lambda - lambda_old;
55: /* check for convergence */
56: if (PetscAbsReal(delLambda) < steptol*lambda) break;
57: if (PetscAbsScalar(fty) / PetscAbsScalar(fty_init) < rtol) break;
58: if (PetscAbsScalar(fty) < atol && i > 0) break;
59: if (monitor) {
60: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
61: PetscViewerASCIIPrintf(monitor," Line search: lambdas = [%g, %g], ftys = [%g, %g]\n",(double)lambda, (double)lambda_old, (double)PetscRealPart(fty), (double)PetscRealPart(fty_old));
62: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
63: }
65: /* compute the search direction */
66: if (linesearch->order == SNES_LINESEARCH_ORDER_LINEAR) {
67: s = (fty - fty_old) / delLambda;
68: } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
69: VecCopy(X, W);
70: VecAXPY(W, -0.5*(lambda + lambda_old), Y);
71: if (linesearch->ops->viproject) {
72: (*linesearch->ops->viproject)(snes, W);
73: }
74: (*linesearch->ops->snesfunc)(snes,W,F);
75: VecDot(F, Y, &fty_mid1);
76: s = (3.*fty - 4.*fty_mid1 + fty_old) / delLambda;
77: } else {
78: VecCopy(X, W);
79: VecAXPY(W, -0.5*(lambda + lambda_old), Y);
80: if (linesearch->ops->viproject) {
81: (*linesearch->ops->viproject)(snes, W);
82: }
83: (*linesearch->ops->snesfunc)(snes,W,F);
84: VecDot(F, Y, &fty_mid1);
85: VecCopy(X, W);
86: VecAXPY(W, -(lambda + 0.5*(lambda - lambda_old)), Y);
87: if (linesearch->ops->viproject) {
88: (*linesearch->ops->viproject)(snes, W);
89: }
90: (*linesearch->ops->snesfunc)(snes, W, F);
91: VecDot(F, Y, &fty_mid2);
92: s = (2.*fty_mid2 + 3.*fty - 6.*fty_mid1 + fty_old) / (3.*delLambda);
93: }
94: /* if the solve is going in the wrong direction, fix it */
95: if (PetscRealPart(s) > 0.) s = -s;
96: if (s == 0.0) break;
97: lambda_update = lambda - PetscRealPart(fty / s);
99: /* switch directions if we stepped out of bounds */
100: if (lambda_update < steptol) lambda_update = lambda + PetscRealPart(fty / s);
102: if (PetscIsInfOrNanReal(lambda_update)) break;
103: if (lambda_update > maxstep) break;
105: /* compute the new state of the line search */
106: lambda_old = lambda;
107: lambda = lambda_update;
108: fty_old = fty;
109: }
110: /* construct the solution */
111: VecCopy(X, W);
112: VecAXPY(W, -lambda, Y);
113: if (linesearch->ops->viproject) {
114: (*linesearch->ops->viproject)(snes, W);
115: }
116: /* postcheck */
117: SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w);
118: if (changed_y) {
119: VecAXPY(X, -lambda, Y);
120: if (linesearch->ops->viproject) {
121: (*linesearch->ops->viproject)(snes, X);
122: }
123: } else {
124: VecCopy(W, X);
125: }
126: (*linesearch->ops->snesfunc)(snes,X,F);
128: SNESLineSearchComputeNorms(linesearch);
129: SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);
131: SNESLineSearchSetLambda(linesearch, lambda);
133: if (monitor) {
134: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
135: PetscViewerASCIIPrintf(monitor," Line search terminated: lambda = %g, fnorms = %g\n", (double)lambda, (double)gnorm);
136: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
137: }
138: if (lambda <= steptol) {
139: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);
140: }
141: return(0);
142: }
144: /*MC
145: SNESLINESEARCHCP - Critical point line search. This line search assumes that there exists some
146: artificial G(x) for which the SNESFunction F(x) = grad G(x). Therefore, this line search seeks
147: to find roots of dot(F, Y) via a secant method.
149: Options Database Keys:
150: + -snes_linesearch_minlambda <minlambda> - the minimum acceptable lambda
151: . -snes_linesearch_maxstep <length> - the algorithm insures that a step length is never longer than this value
152: . -snes_linesearch_damping <damping> - initial trial step length is scaled by this factor, default is 1.0
153: - -snes_linesearch_max_it <max_it> - the maximum number of secant steps performed.
155: Notes:
156: This method does NOT use the objective function if it is provided with SNESSetObjective().
158: This method is the preferred line search for SNESQN and SNESNCG.
160: Level: advanced
162: .keywords: SNES, SNESLineSearch, damping
164: .seealso: SNESLineSearchCreate(), SNESLineSearchSetType()
165: M*/
166: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_CP(SNESLineSearch linesearch)
167: {
169: linesearch->ops->apply = SNESLineSearchApply_CP;
170: linesearch->ops->destroy = NULL;
171: linesearch->ops->setfromoptions = NULL;
172: linesearch->ops->reset = NULL;
173: linesearch->ops->view = NULL;
174: linesearch->ops->setup = NULL;
175: linesearch->order = SNES_LINESEARCH_ORDER_LINEAR;
177: linesearch->max_its = 1;
178: return(0);
179: }