Actual source code: gpcglinesearch.c
1: #include "private/taolinesearch_impl.h"
2: #include "gpcglinesearch.h"
4: /* ---------------------------------------------------------- */
8: static PetscErrorCode TaoLineSearchDestroy_GPCG(TaoLineSearch ls)
9: {
10: PetscErrorCode ierr;
11: TAOLINESEARCH_GPCG_CTX *ctx = (TAOLINESEARCH_GPCG_CTX *)ls->data;
14: if (ctx->W1) {
15: VecDestroy(&ctx->W1);
16: }
17: if (ctx->W2) {
18: VecDestroy(&ctx->W2);
19: }
20: if (ctx->Gold) {
21: VecDestroy(&ctx->Gold);
22: }
23: if (ctx->x) {
24: PetscObjectDereference((PetscObject)ctx->x);
25: }
27: PetscFree(ctx);
28: ls->data = 0;
29: return(0);
30: }
33: /*------------------------------------------------------------*/
36: static PetscErrorCode TaoLineSearchView_GPCG(TaoLineSearch ls, PetscViewer viewer)
37: {
38: PetscBool isascii;
39: PetscErrorCode ierr;
42: PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
44: if (isascii) {
45: PetscViewerASCIIPrintf(viewer," GPCG Line search");
46: } else {
47: SETERRQ1(((PetscObject)ls)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for GPCG LineSearch",((PetscObject)viewer)->type_name);
48: }
49: return(0);
50: }
52: /*------------------------------------------------------------*/
55: static PetscErrorCode TaoLineSearchApply_GPCG(TaoLineSearch ls, Vec x,
56: PetscReal *f, Vec g, Vec s)
58: {
59: TAOLINESEARCH_GPCG_CTX *neP = (TAOLINESEARCH_GPCG_CTX *)ls->data;
60: PetscErrorCode ierr;
61: PetscInt i;
62: PetscBool g_computed=PETSC_FALSE; /* to prevent extra gradient computation */
63: PetscReal d1,finit,actred,prered,rho, gdx;
66: /* ls->stepmin - lower bound for step */
67: /* ls->stepmax - upper bound for step */
68: /* ls->rtol - relative tolerance for an acceptable step */
69: /* ls->ftol - tolerance for sufficient decrease condition */
70: /* ls->gtol - tolerance for curvature condition */
71: /* ls->nfeval - number of function evaluations */
72: /* ls->nfeval - number of function/gradient evaluations */
73: /* ls->max_funcs - maximum number of function evaluations */
75: ls->reason = TAOLINESEARCH_CONTINUE_ITERATING;
76: ls->step = ls->initstep;
77: if (!neP->W2) {
78: VecDuplicate(x,&neP->W2);
79: VecDuplicate(x,&neP->W1);
80: VecDuplicate(x,&neP->Gold);
81: neP->x = x;
82: PetscObjectReference((PetscObject)neP->x);
83: }
85: /* If X has changed, remake work vectors */
86: else if (x != neP->x) {
87: VecDestroy(&neP->x);
88: VecDestroy(&neP->W1);
89: VecDestroy(&neP->W2);
90: VecDestroy(&neP->Gold);
91: VecDuplicate(x,&neP->W1);
92: VecDuplicate(x,&neP->W2);
93: VecDuplicate(x,&neP->Gold);
94: PetscObjectDereference((PetscObject)neP->x);
95: neP->x = x;
96: PetscObjectReference((PetscObject)neP->x);
97: }
99: VecDot(g,s,&gdx);
100: VecCopy(x,neP->W2);
101: VecCopy(g,neP->Gold);
102: if (ls->bounded) {
103: /* Compute the smallest steplength that will make one nonbinding variable
104: equal the bound */
105: VecStepBoundInfo(x,ls->lower,ls->upper,s,&rho,&actred,&d1);
106: ls->step = PetscMin(ls->step,d1);
107: }
108: rho=0; actred=0;
110: if (ls->step < 0) {
111: PetscInfo1(ls,"Line search error: initial step parameter %G < 0\n",ls->step);
112: ls->reason = TAOLINESEARCH_HALTED_OTHER;
113: return(0);
114: }
116: /* Initialization */
117: finit = *f;
118: for (i=0; i< ls->max_funcs; i++) {
119:
120: /* Force the step to be within the bounds */
121: ls->step = PetscMax(ls->step,ls->stepmin);
122: ls->step = PetscMin(ls->step,ls->stepmax);
124: VecCopy(x,neP->W2);
125: VecAXPY(neP->W2,ls->step,s);
126: if (ls->bounded) {
127: VecMedian(ls->lower,neP->W2,ls->upper,neP->W2);
128: }
130: /* Gradient is not needed here. Unless there is a separate
131: gradient routine, compute it here anyway to prevent recomputing at
132: the end of the line search */
133: if (ls->hasobjective) {
134: TaoLineSearchComputeObjective(ls,neP->W2,f);
135: g_computed=PETSC_FALSE;
136: } else if (ls->usegts){
137: TaoLineSearchComputeObjectiveAndGTS(ls,neP->W2,f,&gdx);
138: g_computed=PETSC_FALSE;
139: } else {
140: TaoLineSearchComputeObjectiveAndGradient(ls,neP->W2,f,g);
141: g_computed=PETSC_TRUE;
142: }
144: if (0 == i) {
145: ls->f_fullstep = *f;
146: }
148: actred = *f - finit;
149: VecCopy(neP->W2,neP->W1);
150: VecAXPY(neP->W1,-1.0,x); /* W1 = W2 - X */
151: VecDot(neP->W1,neP->Gold,&prered);
152:
153: if (fabs(prered)<1.0e-100) prered=1.0e-12;
154: rho = actred/prered;
155: /*
156: If sufficient progress has been obtained, accept the
157: point. Otherwise, backtrack.
158: */
160: if (rho > ls->ftol){
161: break;
162: } else{
163: ls->step = (ls->step)/2;
164: }
166: /* Convergence testing */
167:
168: if (ls->step <= ls->stepmin || ls->step >= ls->stepmax) {
169: ls->reason = TAOLINESEARCH_HALTED_OTHER;
170: PetscInfo(ls,"Rounding errors may prevent further progress. May not be a step satisfying\n");
171: PetscInfo(ls,"sufficient decrease and curvature conditions. Tolerances may be too small.\n");
172: break;
173: }
174: if (ls->step == ls->stepmax) {
175: PetscInfo1(ls,"Step is at the upper bound, stepmax (%G)\n",ls->stepmax);
176: ls->reason = TAOLINESEARCH_HALTED_UPPERBOUND;
177: break;
178: }
179: if (ls->step == ls->stepmin) {
180: PetscInfo1(ls,"Step is at the lower bound, stepmin (%G)\n",ls->stepmin);
181: ls->reason = TAOLINESEARCH_HALTED_LOWERBOUND;
182: break;
183: }
184: if ((ls->nfeval+ls->nfgeval) >= ls->max_funcs) {
185: PetscInfo2(ls,"Number of line search function evals (%D) > maximum (%D)\n",ls->nfeval+ls->nfgeval,ls->max_funcs);
186: ls->reason = TAOLINESEARCH_HALTED_MAXFCN;
187: break;
188: }
189: if ((neP->bracket) && (ls->stepmax - ls->stepmin <= ls->rtol*ls->stepmax)){
190: PetscInfo1(ls,"Relative width of interval of uncertainty is at most rtol (%G)\n",ls->rtol);
191: ls->reason = TAOLINESEARCH_HALTED_RTOL;
192: break;
193: }
194: }
195: PetscInfo2(ls,"%D function evals in line search, step = %G\n",ls->nfeval+ls->nfgeval,ls->step);
196: /* set new solution vector and compute gradient if necessary */
197: VecCopy(neP->W2, x);
198: if (!g_computed) {
199: TaoLineSearchComputeGradient(ls,x,g);
200: }
201: return(0);
202: }
204: /* ---------------------------------------------------------- */
208: PetscErrorCode TaoLineSearchCreate_GPCG(TaoLineSearch ls)
209: {
211: TAOLINESEARCH_GPCG_CTX *neP;
217: ls->ftol = 0.05;
218: ls->rtol = 0.0;
219: ls->gtol = 0.0;
220: ls->stepmin = 1.0e-20;
221: ls->stepmax = 1.0e+20;
222: ls->nfeval = 0;
223: ls->max_funcs = 30;
224: ls->step = 1.0;
226: PetscNewLog(ls,TAOLINESEARCH_GPCG_CTX,&neP);
227: neP->bracket = 0;
228: neP->infoc = 1;
229: ls->data = (void*)neP;
231: ls->ops->setup = 0;
232: ls->ops->apply=TaoLineSearchApply_GPCG;
233: ls->ops->view =TaoLineSearchView_GPCG;
234: ls->ops->destroy=TaoLineSearchDestroy_GPCG;
235: ls->ops->setfromoptions=0;
237: return(0);
238: }