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