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