xref: /petsc/src/snes/linesearch/impls/cp/linesearchcp.c (revision ea5d4fccf296dd2bbd0f9c3a3343651cb1066da7)
1 #include <private/linesearchimpl.h>
2 #include <petscsnes.h>
3 
4 #undef __FUNCT__
5 #define __FUNCT__ "PetscLineSearchApply_CP"
6 
7 /*@C
8    LineSearchCP - This routine is not a line search at all;
9    it simply uses the full step.  Thus, this routine is intended
10    to serve as a template and is not recommended for general use.
11 
12    Logically Collective on SNES and Vec
13 
14    Input Parameters:
15 +  snes - nonlinear context
16 .  lsctx - optional context for line search (not used here)
17 .  x - current iterate
18 .  f - residual evaluated at x
19 .  y - search direction
20 .  fnorm - 2-norm of f
21 -  xnorm - norm of x if known, otherwise 0
22 
23    Output Parameters:
24 +  g - residual evaluated at new iterate y
25 .  w - new iterate
26 .  gnorm - 2-norm of g
27 .  ynorm - 2-norm of search length
28 -  flag - PETSC_TRUE on success, PETSC_FALSE on failure
29 
30    Options Database Key:
31 .  -snes_ls basic - Activates SNESLineSearchNo()
32 
33    Level: advanced
34 
35 .keywords: SNES, nonlinear, line search, cubic
36 
37 .seealso: SNESLineSearchCubic(), SNESLineSearchQuadratic(),
38           SNESLineSearchSet(), SNESLineSearchNoNorms()
39 @*/
40 PetscErrorCode  PetscLineSearchApply_CP(PetscLineSearch linesearch)
41 {
42 
43   PetscBool      changed_y, changed_w, domainerror;
44   PetscErrorCode ierr;
45   Vec             X, Y, F, W;
46   SNES            snes;
47   PetscReal       xnorm, ynorm, gnorm, steptol, atol, rtol, ltol, maxstep;
48 
49   PetscReal       lambda, lambda_old, lambda_update, delLambda;
50   PetscScalar     fty, fty_old;
51   PetscInt        i, max_its;
52 
53   PetscViewer     monitor;
54 
55   PetscFunctionBegin;
56 
57   ierr = PetscLineSearchGetVecs(linesearch, &X, &F, &Y, &W, PETSC_NULL);CHKERRQ(ierr);
58   ierr = PetscLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);CHKERRQ(ierr);
59   ierr = PetscLineSearchGetSNES(linesearch, &snes);CHKERRQ(ierr);
60   ierr = PetscLineSearchGetLambda(linesearch, &lambda);CHKERRQ(ierr);
61   ierr = PetscLineSearchGetTolerances(linesearch, &steptol, &maxstep, &rtol, &atol, &ltol, &max_its);CHKERRQ(ierr);
62   ierr = PetscLineSearchSetSuccess(linesearch, PETSC_TRUE);CHKERRQ(ierr);
63   ierr = PetscLineSearchGetMonitor(linesearch, &monitor);CHKERRQ(ierr);
64 
65   /* precheck */
66   ierr = PetscLineSearchPreCheck(linesearch, &changed_y);CHKERRQ(ierr);
67   lambda_old = 0.0;
68   ierr = VecDot(F, Y, &fty_old);CHKERRQ(ierr);
69   for (i = 0; i < max_its; i++) {
70 
71     /* compute the norm at lambda */
72     ierr = VecCopy(X, W);CHKERRQ(ierr);
73     ierr = VecAXPY(W, -lambda, Y);CHKERRQ(ierr);
74     if (linesearch->ops->viproject) {
75       ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr);
76     }
77     ierr = SNESComputeFunction(snes, W, F);CHKERRQ(ierr);
78 
79     ierr = VecDot(F, Y, &fty);CHKERRQ(ierr);
80 
81     delLambda    = lambda - lambda_old;
82     if (PetscAbsReal(delLambda) < steptol) break;
83     if (monitor) {
84       ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
85       ierr = PetscViewerASCIIPrintf(monitor,"    Line search: lambdas = [%g, %g], ftys = [%g, %g]\n",
86                                     lambda, lambda_old, PetscRealPart(fty), PetscRealPart(fty_old));CHKERRQ(ierr);
87       ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
88     }
89 
90     /* compute the search direction */
91     lambda_update =  PetscRealPart((fty*lambda_old - fty_old*lambda) / (fty - fty_old));
92     if (PetscIsInfOrNanScalar(lambda_update)) break;
93     if (lambda_update > maxstep) {
94       break;
95     }
96 
97     /* compute the new state of the line search */
98     lambda_old = lambda;
99     lambda = lambda_update;
100     fty_old = fty;
101   }
102   /* construct the solution */
103   ierr = VecCopy(X, W);CHKERRQ(ierr);
104   ierr = VecAXPY(W, -lambda, Y);CHKERRQ(ierr);
105   if (linesearch->ops->viproject) {
106     ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr);
107   }
108   /* postcheck */
109   ierr = PetscLineSearchPostCheck(linesearch, &changed_y, &changed_w);CHKERRQ(ierr);
110   if (changed_y) {
111     ierr = VecAXPY(X, -lambda, Y);CHKERRQ(ierr);
112     if (linesearch->ops->viproject) {
113       ierr = (*linesearch->ops->viproject)(snes, X);CHKERRQ(ierr);
114     }
115   } else {
116     ierr = VecCopy(W, X);CHKERRQ(ierr);
117   }
118   ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr);
119   ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr);
120   if (domainerror) {
121     ierr = PetscLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
122     PetscFunctionReturn(0);
123   }
124 
125   ierr = PetscLineSearchComputeNorms(linesearch);CHKERRQ(ierr);
126   ierr = PetscLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);CHKERRQ(ierr);
127 
128   ierr = PetscLineSearchSetLambda(linesearch, lambda);CHKERRQ(ierr);
129 
130   if (monitor) {
131     ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
132     ierr = PetscViewerASCIIPrintf(monitor,"    Line search terminated: lambda = %g, fnorms = %g\n", lambda, gnorm);CHKERRQ(ierr);
133     ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
134   }
135   if (lambda <= steptol) {
136     ierr = PetscLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
137   }
138   PetscFunctionReturn(0);
139 }
140 
141 EXTERN_C_BEGIN
142 #undef __FUNCT__
143 #define __FUNCT__ "PetscLineSearchCreate_CP"
144 PetscErrorCode PetscLineSearchCreate_CP(PetscLineSearch linesearch)
145 {
146   PetscFunctionBegin;
147   linesearch->ops->apply          = PetscLineSearchApply_CP;
148   linesearch->ops->destroy        = PETSC_NULL;
149   linesearch->ops->setfromoptions = PETSC_NULL;
150   linesearch->ops->reset          = PETSC_NULL;
151   linesearch->ops->view           = PETSC_NULL;
152   linesearch->ops->setup          = PETSC_NULL;
153   PetscFunctionReturn(0);
154 }
155 EXTERN_C_END
156