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: }