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