Actual source code: gpcglinesearch.c
1: #include <petsc/private/taolinesearchimpl.h>
2: #include <../src/tao/linesearch/impls/gpcglinesearch/gpcglinesearch.h>
4: /* ---------------------------------------------------------- */
6: static PetscErrorCode TaoLineSearchDestroy_GPCG(TaoLineSearch ls)
7: {
8: TaoLineSearch_GPCG *ctx = (TaoLineSearch_GPCG *)ls->data;
10: VecDestroy(&ctx->W1);
11: VecDestroy(&ctx->W2);
12: VecDestroy(&ctx->Gold);
13: VecDestroy(&ctx->x);
14: PetscFree(ls->data);
15: return 0;
16: }
18: /*------------------------------------------------------------*/
19: static PetscErrorCode TaoLineSearchView_GPCG(TaoLineSearch ls, PetscViewer viewer)
20: {
21: PetscBool isascii;
23: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
24: if (isascii) {
25: PetscViewerASCIIPrintf(viewer," GPCG Line search");
26: }
27: return 0;
28: }
30: /*------------------------------------------------------------*/
31: static PetscErrorCode TaoLineSearchApply_GPCG(TaoLineSearch ls, Vec x, PetscReal *f, Vec g, Vec s)
32: {
33: TaoLineSearch_GPCG *neP = (TaoLineSearch_GPCG *)ls->data;
34: PetscInt i;
35: PetscBool g_computed=PETSC_FALSE; /* to prevent extra gradient computation */
36: PetscReal d1,finit,actred,prered,rho, gdx;
38: /* ls->stepmin - lower bound for step */
39: /* ls->stepmax - upper bound for step */
40: /* ls->rtol - relative tolerance for an acceptable step */
41: /* ls->ftol - tolerance for sufficient decrease condition */
42: /* ls->gtol - tolerance for curvature condition */
43: /* ls->nfeval - number of function evaluations */
44: /* ls->nfeval - number of function/gradient evaluations */
45: /* ls->max_funcs - maximum number of function evaluations */
47: TaoLineSearchMonitor(ls, 0, *f, 0.0);
49: ls->reason = TAOLINESEARCH_CONTINUE_ITERATING;
50: ls->step = ls->initstep;
51: if (!neP->W2) {
52: VecDuplicate(x,&neP->W2);
53: VecDuplicate(x,&neP->W1);
54: VecDuplicate(x,&neP->Gold);
55: neP->x = x;
56: PetscObjectReference((PetscObject)neP->x);
57: } else if (x != neP->x) {
58: VecDestroy(&neP->x);
59: VecDestroy(&neP->W1);
60: VecDestroy(&neP->W2);
61: VecDestroy(&neP->Gold);
62: VecDuplicate(x,&neP->W1);
63: VecDuplicate(x,&neP->W2);
64: VecDuplicate(x,&neP->Gold);
65: PetscObjectDereference((PetscObject)neP->x);
66: neP->x = x;
67: PetscObjectReference((PetscObject)neP->x);
68: }
70: VecDot(g,s,&gdx);
71: if (gdx > 0) {
72: PetscInfo(ls,"Line search error: search direction is not descent direction. dot(g,s) = %g\n",(double)gdx);
73: ls->reason = TAOLINESEARCH_FAILED_ASCENT;
74: return 0;
75: }
76: VecCopy(x,neP->W2);
77: VecCopy(g,neP->Gold);
78: if (ls->bounded) {
79: /* Compute the smallest steplength that will make one nonbinding variable equal the bound */
80: VecStepBoundInfo(x,s,ls->lower,ls->upper,&rho,&actred,&d1);
81: ls->step = PetscMin(ls->step,d1);
82: }
83: rho=0; actred=0;
85: if (ls->step < 0) {
86: PetscInfo(ls,"Line search error: initial step parameter %g< 0\n",(double)ls->step);
87: ls->reason = TAOLINESEARCH_HALTED_OTHER;
88: return 0;
89: }
91: /* Initialization */
92: finit = *f;
93: for (i=0; i< ls->max_funcs; i++) {
94: /* Force the step to be within the bounds */
95: ls->step = PetscMax(ls->step,ls->stepmin);
96: ls->step = PetscMin(ls->step,ls->stepmax);
98: VecCopy(x,neP->W2);
99: VecAXPY(neP->W2,ls->step,s);
100: if (ls->bounded) {
101: /* Make sure new vector is numerically within bounds */
102: VecMedian(neP->W2,ls->lower,ls->upper,neP->W2);
103: }
105: /* Gradient is not needed here. Unless there is a separate
106: gradient routine, compute it here anyway to prevent recomputing at
107: the end of the line search */
108: if (ls->hasobjective) {
109: TaoLineSearchComputeObjective(ls,neP->W2,f);
110: g_computed=PETSC_FALSE;
111: } else if (ls->usegts) {
112: TaoLineSearchComputeObjectiveAndGTS(ls,neP->W2,f,&gdx);
113: g_computed=PETSC_FALSE;
114: } else {
115: TaoLineSearchComputeObjectiveAndGradient(ls,neP->W2,f,g);
116: g_computed=PETSC_TRUE;
117: }
119: TaoLineSearchMonitor(ls, i+1, *f, ls->step);
121: if (0 == i) {
122: ls->f_fullstep = *f;
123: }
125: actred = *f - finit;
126: VecCopy(neP->W2,neP->W1);
127: VecAXPY(neP->W1,-1.0,x); /* W1 = W2 - X */
128: VecDot(neP->W1,neP->Gold,&prered);
130: if (PetscAbsReal(prered)<1.0e-100) prered=1.0e-12;
131: rho = actred/prered;
133: /*
134: If sufficient progress has been obtained, accept the
135: point. Otherwise, backtrack.
136: */
138: if (actred > 0) {
139: PetscInfo(ls,"Step resulted in ascent, rejecting.\n");
140: ls->step = (ls->step)/2;
141: } else if (rho > ls->ftol) {
142: break;
143: } else{
144: ls->step = (ls->step)/2;
145: }
147: /* Convergence testing */
149: if (ls->step <= ls->stepmin || ls->step >= ls->stepmax) {
150: ls->reason = TAOLINESEARCH_HALTED_OTHER;
151: PetscInfo(ls,"Rounding errors may prevent further progress. May not be a step satisfying\n");
152: PetscInfo(ls,"sufficient decrease and curvature conditions. Tolerances may be too small.\n");
153: break;
154: }
155: if (ls->step == ls->stepmax) {
156: PetscInfo(ls,"Step is at the upper bound, stepmax (%g)\n",(double)ls->stepmax);
157: ls->reason = TAOLINESEARCH_HALTED_UPPERBOUND;
158: break;
159: }
160: if (ls->step == ls->stepmin) {
161: PetscInfo(ls,"Step is at the lower bound, stepmin (%g)\n",(double)ls->stepmin);
162: ls->reason = TAOLINESEARCH_HALTED_LOWERBOUND;
163: break;
164: }
165: if ((ls->nfeval+ls->nfgeval) >= ls->max_funcs) {
166: PetscInfo(ls,"Number of line search function evals (%D) > maximum (%D)\n",ls->nfeval+ls->nfgeval,ls->max_funcs);
167: ls->reason = TAOLINESEARCH_HALTED_MAXFCN;
168: break;
169: }
170: if ((neP->bracket) && (ls->stepmax - ls->stepmin <= ls->rtol*ls->stepmax)) {
171: PetscInfo(ls,"Relative width of interval of uncertainty is at most rtol (%g)\n",(double)ls->rtol);
172: ls->reason = TAOLINESEARCH_HALTED_RTOL;
173: break;
174: }
175: }
176: PetscInfo(ls,"%D function evals in line search, step = %g\n",ls->nfeval+ls->nfgeval,(double)ls->step);
177: /* set new solution vector and compute gradient if necessary */
178: VecCopy(neP->W2, x);
179: if (ls->reason == TAOLINESEARCH_CONTINUE_ITERATING) {
180: ls->reason = TAOLINESEARCH_SUCCESS;
181: }
182: if (!g_computed) {
183: TaoLineSearchComputeGradient(ls,x,g);
184: }
185: return 0;
186: }
188: /* ---------------------------------------------------------- */
190: /*MC
191: TAOLINESEARCHGPCG - Special line-search method for the Gradient-Projected Conjugate Gradient (TAOGPCG) algorithm.
192: Should not be used with any other algorithm.
194: Level: developer
196: .keywords: Tao, linesearch
197: M*/
198: PETSC_EXTERN PetscErrorCode TaoLineSearchCreate_GPCG(TaoLineSearch ls)
199: {
200: TaoLineSearch_GPCG *neP;
202: ls->ftol = 0.05;
203: ls->rtol = 0.0;
204: ls->gtol = 0.0;
205: ls->stepmin = 1.0e-20;
206: ls->stepmax = 1.0e+20;
207: ls->nfeval = 0;
208: ls->max_funcs = 30;
209: ls->step = 1.0;
211: PetscNewLog(ls,&neP);
212: neP->bracket = 0;
213: neP->infoc = 1;
214: ls->data = (void*)neP;
216: ls->ops->setup = NULL;
217: ls->ops->reset = NULL;
218: ls->ops->apply = TaoLineSearchApply_GPCG;
219: ls->ops->view = TaoLineSearchView_GPCG;
220: ls->ops->destroy = TaoLineSearchDestroy_GPCG;
221: ls->ops->setfromoptions = NULL;
222: ls->ops->monitor = NULL;
223: return 0;
224: }