xref: /petsc/src/tao/unconstrained/impls/cg/taocg.c (revision e4cb33bb7dbdbae9285fba102465ca0f1dcb3977) 
1 #include <petsctaolinesearch.h>
2 #include <../src/tao/unconstrained/impls/cg/taocg.h>
3 
4 #define CG_FletcherReeves       0
5 #define CG_PolakRibiere         1
6 #define CG_PolakRibierePlus     2
7 #define CG_HestenesStiefel      3
8 #define CG_DaiYuan              4
9 #define CG_Types                5
10 
11 static const char *CG_Table[64] = {"fr", "pr", "prp", "hs", "dy"};
12 
13 #undef __FUNCT__
14 #define __FUNCT__ "TaoSolve_CG"
15 static PetscErrorCode TaoSolve_CG(Tao tao)
16 {
17   TAO_CG                       *cgP = (TAO_CG*)tao->data;
18   PetscErrorCode               ierr;
19   TaoConvergedReason           reason = TAO_CONTINUE_ITERATING;
20   TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;
21   PetscReal                    step=1.0,f,gnorm,gnorm2,delta,gd,ginner,beta;
22   PetscReal                    gd_old,gnorm2_old,f_old;
23   PetscInt                     iter=0;
24 
25   PetscFunctionBegin;
26   if (tao->XL || tao->XU || tao->ops->computebounds) {
27     ierr = PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by cg algorithm\n");CHKERRQ(ierr);
28   }
29 
30   /*  Check convergence criteria */
31   ierr = TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);CHKERRQ(ierr);
32   ierr = VecNorm(tao->gradient,NORM_2,&gnorm);CHKERRQ(ierr);
33   if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
34 
35   ierr = TaoMonitor(tao, iter, f, gnorm, 0.0, step, &reason);CHKERRQ(ierr);
36   if (reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0);
37 
38   /*  Set initial direction to -gradient */
39   ierr = VecCopy(tao->gradient, tao->stepdirection);CHKERRQ(ierr);
40   ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr);
41   gnorm2 = gnorm*gnorm;
42 
43   /*  Set initial scaling for the function */
44   if (f != 0.0) {
45     delta = 2.0*PetscAbsScalar(f) / gnorm2;
46     delta = PetscMax(delta,cgP->delta_min);
47     delta = PetscMin(delta,cgP->delta_max);
48   } else {
49     delta = 2.0 / gnorm2;
50     delta = PetscMax(delta,cgP->delta_min);
51     delta = PetscMin(delta,cgP->delta_max);
52   }
53   /*  Set counter for gradient and reset steps */
54   cgP->ngradsteps = 0;
55   cgP->nresetsteps = 0;
56 
57   while (1) {
58     /*  Save the current gradient information */
59     f_old = f;
60     gnorm2_old = gnorm2;
61     ierr = VecCopy(tao->solution, cgP->X_old);CHKERRQ(ierr);
62     ierr = VecCopy(tao->gradient, cgP->G_old);CHKERRQ(ierr);
63     ierr = VecDot(tao->gradient, tao->stepdirection, &gd);CHKERRQ(ierr);
64     if ((gd >= 0) || PetscIsInfOrNanReal(gd)) {
65       ++cgP->ngradsteps;
66       if (f != 0.0) {
67         delta = 2.0*PetscAbsScalar(f) / gnorm2;
68         delta = PetscMax(delta,cgP->delta_min);
69         delta = PetscMin(delta,cgP->delta_max);
70       } else {
71         delta = 2.0 / gnorm2;
72         delta = PetscMax(delta,cgP->delta_min);
73         delta = PetscMin(delta,cgP->delta_max);
74       }
75 
76       ierr = VecCopy(tao->gradient, tao->stepdirection);CHKERRQ(ierr);
77       ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr);
78     }
79 
80     /*  Search direction for improving point */
81     ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,delta);
82     ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status);CHKERRQ(ierr);
83     ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr);
84     if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
85       /*  Linesearch failed */
86       /*  Reset factors and use scaled gradient step */
87       ++cgP->nresetsteps;
88       f = f_old;
89       gnorm2 = gnorm2_old;
90       ierr = VecCopy(cgP->X_old, tao->solution);CHKERRQ(ierr);
91       ierr = VecCopy(cgP->G_old, tao->gradient);CHKERRQ(ierr);
92 
93       if (f != 0.0) {
94         delta = 2.0*PetscAbsScalar(f) / gnorm2;
95         delta = PetscMax(delta,cgP->delta_min);
96         delta = PetscMin(delta,cgP->delta_max);
97       } else {
98         delta = 2.0 / gnorm2;
99         delta = PetscMax(delta,cgP->delta_min);
100         delta = PetscMin(delta,cgP->delta_max);
101       }
102 
103       ierr = VecCopy(tao->gradient, tao->stepdirection);CHKERRQ(ierr);
104       ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr);
105 
106       ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,delta);
107       ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status);CHKERRQ(ierr);
108       ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr);
109 
110       if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
111         /*  Linesearch failed again */
112         /*  switch to unscaled gradient */
113         f = f_old;
114         gnorm2 = gnorm2_old;
115         ierr = VecCopy(cgP->X_old, tao->solution);CHKERRQ(ierr);
116         ierr = VecCopy(cgP->G_old, tao->gradient);CHKERRQ(ierr);
117         delta = 1.0;
118         ierr = VecCopy(tao->solution, tao->stepdirection);CHKERRQ(ierr);
119         ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr);
120 
121         ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,delta);
122         ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status);CHKERRQ(ierr);
123         ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr);
124         if (ls_status != TAOLINESEARCH_SUCCESS &&
125             ls_status != TAOLINESEARCH_SUCCESS_USER) {
126 
127           /*  Line search failed for last time -- give up */
128           f = f_old;
129           gnorm2 = gnorm2_old;
130           ierr = VecCopy(cgP->X_old, tao->solution);CHKERRQ(ierr);
131           ierr = VecCopy(cgP->G_old, tao->gradient);CHKERRQ(ierr);
132           step = 0.0;
133           reason = TAO_DIVERGED_LS_FAILURE;
134           tao->reason = TAO_DIVERGED_LS_FAILURE;
135         }
136       }
137     }
138 
139     /*  Check for bad value */
140     ierr = VecNorm(tao->gradient,NORM_2,&gnorm);CHKERRQ(ierr);
141     if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1,"User-provided compute function generated Inf or NaN");
142 
143     /*  Check for termination */
144     gnorm2 =gnorm * gnorm;
145     iter++;
146     ierr = TaoMonitor(tao, iter, f, gnorm, 0.0, step, &reason);CHKERRQ(ierr);
147     if (reason != TAO_CONTINUE_ITERATING) {
148       break;
149     }
150 
151     /*  Check for restart condition */
152     ierr = VecDot(tao->gradient, cgP->G_old, &ginner);CHKERRQ(ierr);
153     if (PetscAbsScalar(ginner) >= cgP->eta * gnorm2) {
154       /*  Gradients far from orthognal; use steepest descent direction */
155       beta = 0.0;
156     } else {
157       /*  Gradients close to orthogonal; use conjugate gradient formula */
158       switch (cgP->cg_type) {
159       case CG_FletcherReeves:
160         beta = gnorm2 / gnorm2_old;
161         break;
162 
163       case CG_PolakRibiere:
164         beta = (gnorm2 - ginner) / gnorm2_old;
165         break;
166 
167       case CG_PolakRibierePlus:
168         beta = PetscMax((gnorm2-ginner)/gnorm2_old, 0.0);
169         break;
170 
171       case CG_HestenesStiefel:
172         ierr = VecDot(tao->gradient, tao->stepdirection, &gd);CHKERRQ(ierr);
173         ierr = VecDot(cgP->G_old, tao->stepdirection, &gd_old);CHKERRQ(ierr);
174         beta = (gnorm2 - ginner) / (gd - gd_old);
175         break;
176 
177       case CG_DaiYuan:
178         ierr = VecDot(tao->gradient, tao->stepdirection, &gd);CHKERRQ(ierr);
179         ierr = VecDot(cgP->G_old, tao->stepdirection, &gd_old);CHKERRQ(ierr);
180         beta = gnorm2 / (gd - gd_old);
181         break;
182 
183       default:
184         beta = 0.0;
185         break;
186       }
187     }
188 
189     /*  Compute the direction d=-g + beta*d */
190     ierr = VecAXPBY(tao->stepdirection, -1.0, beta, tao->gradient);CHKERRQ(ierr);
191 
192     /*  update initial steplength choice */
193     delta = 1.0;
194     delta = PetscMax(delta, cgP->delta_min);
195     delta = PetscMin(delta, cgP->delta_max);
196   }
197   PetscFunctionReturn(0);
198 }
199 
200 #undef __FUNCT__
201 #define __FUNCT__ "TaoSetUp_CG"
202 static PetscErrorCode TaoSetUp_CG(Tao tao)
203 {
204   TAO_CG         *cgP = (TAO_CG*)tao->data;
205   PetscErrorCode ierr;
206 
207   PetscFunctionBegin;
208   if (!tao->gradient) {ierr = VecDuplicate(tao->solution,&tao->gradient);CHKERRQ(ierr);}
209   if (!tao->stepdirection) {ierr = VecDuplicate(tao->solution,&tao->stepdirection);CHKERRQ(ierr); }
210   if (!cgP->X_old) {ierr = VecDuplicate(tao->solution,&cgP->X_old);CHKERRQ(ierr);}
211   if (!cgP->G_old) {ierr = VecDuplicate(tao->gradient,&cgP->G_old);CHKERRQ(ierr); }
212    PetscFunctionReturn(0);
213 }
214 
215 #undef __FUNCT__
216 #define __FUNCT__ "TaoDestroy_CG"
217 static PetscErrorCode TaoDestroy_CG(Tao tao)
218 {
219   TAO_CG         *cgP = (TAO_CG*) tao->data;
220   PetscErrorCode ierr;
221 
222   PetscFunctionBegin;
223   if (tao->setupcalled) {
224     ierr = VecDestroy(&cgP->X_old);CHKERRQ(ierr);
225     ierr = VecDestroy(&cgP->G_old);CHKERRQ(ierr);
226   }
227   ierr = TaoLineSearchDestroy(&tao->linesearch);CHKERRQ(ierr);
228   ierr = PetscFree(tao->data);CHKERRQ(ierr);
229   PetscFunctionReturn(0);
230 }
231 
232 #undef __FUNCT__
233 #define __FUNCT__ "TaoSetFromOptions_CG"
234 static PetscErrorCode TaoSetFromOptions_CG(Tao tao)
235 {
236    TAO_CG         *cgP = (TAO_CG*)tao->data;
237    PetscErrorCode ierr;
238 
239    PetscFunctionBegin;
240    ierr = TaoLineSearchSetFromOptions(tao->linesearch);CHKERRQ(ierr);
241    ierr = PetscOptionsHead("Nonlinear Conjugate Gradient method for unconstrained optimization");CHKERRQ(ierr);
242    ierr = PetscOptionsReal("-tao_cg_eta","restart tolerance", "", cgP->eta,&cgP->eta,0);CHKERRQ(ierr);
243    ierr = PetscOptionsEList("-tao_cg_type","cg formula", "", CG_Table, CG_Types, CG_Table[cgP->cg_type], &cgP->cg_type, 0);CHKERRQ(ierr);
244    ierr = PetscOptionsReal("-tao_cg_delta_min","minimum delta value", "", cgP->delta_min,&cgP->delta_min,0);CHKERRQ(ierr);
245    ierr = PetscOptionsReal("-tao_cg_delta_max","maximum delta value", "", cgP->delta_max,&cgP->delta_max,0);CHKERRQ(ierr);
246    ierr = PetscOptionsTail();CHKERRQ(ierr);
247    PetscFunctionReturn(0);
248 }
249 
250 #undef __FUNCT__
251 #define __FUNCT__ "TaoView_CG"
252 static PetscErrorCode TaoView_CG(Tao tao, PetscViewer viewer)
253 {
254   PetscBool      isascii;
255   TAO_CG         *cgP = (TAO_CG*)tao->data;
256   PetscErrorCode ierr;
257 
258   PetscFunctionBegin;
259   ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);CHKERRQ(ierr);
260   if (isascii) {
261     ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
262     ierr = PetscViewerASCIIPrintf(viewer, "CG Type: %s\n", CG_Table[cgP->cg_type]);CHKERRQ(ierr);
263     ierr = PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", cgP->ngradsteps);CHKERRQ(ierr);
264     ierr= PetscViewerASCIIPrintf(viewer, "Reset steps: %D\n", cgP->nresetsteps);CHKERRQ(ierr);
265     ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
266   }
267   PetscFunctionReturn(0);
268 }
269 
270 
271 EXTERN_C_BEGIN
272 #undef __FUNCT__
273 #define __FUNCT__ "TaoCreate_CG"
274 PetscErrorCode TaoCreate_CG(Tao tao)
275 {
276   TAO_CG         *cgP;
277   const char     *morethuente_type = TAOLINESEARCH_MT;
278   PetscErrorCode ierr;
279 
280   PetscFunctionBegin;
281   tao->ops->setup = TaoSetUp_CG;
282   tao->ops->solve = TaoSolve_CG;
283   tao->ops->view = TaoView_CG;
284   tao->ops->setfromoptions = TaoSetFromOptions_CG;
285   tao->ops->destroy = TaoDestroy_CG;
286 
287   tao->max_it = 2000;
288   tao->max_funcs = 4000;
289   tao->fatol = 1e-4;
290   tao->frtol = 1e-4;
291 
292   /*  Note: nondefault values should be used for nonlinear conjugate gradient  */
293   /*  method.  In particular, gtol should be less that 0.5; the value used in  */
294   /*  Nocedal and Wright is 0.10.  We use the default values for the  */
295   /*  linesearch because it seems to work better. */
296   ierr = TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);CHKERRQ(ierr);
297   ierr = TaoLineSearchSetType(tao->linesearch, morethuente_type);CHKERRQ(ierr);
298   ierr = TaoLineSearchUseTaoRoutines(tao->linesearch, tao);CHKERRQ(ierr);
299 
300   ierr = PetscNewLog(tao,&cgP);CHKERRQ(ierr);
301   tao->data = (void*)cgP;
302   cgP->eta = 0.1;
303   cgP->delta_min = 1e-7;
304   cgP->delta_max = 100;
305   cgP->cg_type = CG_PolakRibierePlus;
306   PetscFunctionReturn(0);
307 }
308 EXTERN_C_END
309