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