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