Actual source code: blmvm.c

petsc-3.6.1 2015-08-06
Report Typos and Errors
  1: #include <petsctaolinesearch.h>
  2: #include <../src/tao/matrix/lmvmmat.h>
  3: #include <../src/tao/bound/impls/blmvm/blmvm.h>

  5: /*------------------------------------------------------------*/
  8: static PetscErrorCode TaoSolve_BLMVM(Tao tao)
  9: {
 10:   PetscErrorCode               ierr;
 11:   TAO_BLMVM                    *blmP = (TAO_BLMVM *)tao->data;
 12:   TaoConvergedReason           reason = TAO_CONTINUE_ITERATING;
 13:   TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;
 14:   PetscReal                    f, fold, gdx, gnorm;
 15:   PetscReal                    stepsize = 1.0,delta;

 18:   /*  Project initial point onto bounds */
 19:   TaoComputeVariableBounds(tao);
 20:   VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);
 21:   TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);

 23:   /* Check convergence criteria */
 24:   TaoComputeObjectiveAndGradient(tao, tao->solution,&f,blmP->unprojected_gradient);
 25:   VecBoundGradientProjection(blmP->unprojected_gradient,tao->solution, tao->XL,tao->XU,tao->gradient);

 27:   VecNorm(tao->gradient,NORM_2,&gnorm);
 28:   if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf pr NaN");

 30:   TaoMonitor(tao, tao->niter, f, gnorm, 0.0, stepsize, &reason);
 31:   if (reason != TAO_CONTINUE_ITERATING) return(0);

 33:   /* Set initial scaling for the function */
 34:   if (f != 0.0) {
 35:     delta = 2.0*PetscAbsScalar(f) / (gnorm*gnorm);
 36:   } else {
 37:     delta = 2.0 / (gnorm*gnorm);
 38:   }
 39:   MatLMVMSetDelta(blmP->M,delta);

 41:   /* Set counter for gradient/reset steps */
 42:   blmP->grad = 0;
 43:   blmP->reset = 0;

 45:   /* Have not converged; continue with Newton method */
 46:   while (reason == TAO_CONTINUE_ITERATING) {
 47:     /* Compute direction */
 48:     MatLMVMUpdate(blmP->M, tao->solution, tao->gradient);
 49:     MatLMVMSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection);
 50:     VecBoundGradientProjection(tao->stepdirection,tao->solution,tao->XL,tao->XU,tao->gradient);

 52:     /* Check for success (descent direction) */
 53:     VecDot(blmP->unprojected_gradient, tao->gradient, &gdx);
 54:     if (gdx <= 0) {
 55:       /* Step is not descent or solve was not successful
 56:          Use steepest descent direction (scaled) */
 57:       ++blmP->grad;

 59:       if (f != 0.0) {
 60:         delta = 2.0*PetscAbsScalar(f) / (gnorm*gnorm);
 61:       } else {
 62:         delta = 2.0 / (gnorm*gnorm);
 63:       }
 64:       MatLMVMSetDelta(blmP->M,delta);
 65:       MatLMVMReset(blmP->M);
 66:       MatLMVMUpdate(blmP->M, tao->solution, blmP->unprojected_gradient);
 67:       MatLMVMSolve(blmP->M,blmP->unprojected_gradient, tao->stepdirection);
 68:     }
 69:     VecScale(tao->stepdirection,-1.0);

 71:     /* Perform the linesearch */
 72:     fold = f;
 73:     VecCopy(tao->solution, blmP->Xold);
 74:     VecCopy(blmP->unprojected_gradient, blmP->Gold);
 75:     TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
 76:     TaoLineSearchApply(tao->linesearch, tao->solution, &f, blmP->unprojected_gradient, tao->stepdirection, &stepsize, &ls_status);
 77:     TaoAddLineSearchCounts(tao);

 79:     if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
 80:       /* Linesearch failed
 81:          Reset factors and use scaled (projected) gradient step */
 82:       ++blmP->reset;

 84:       f = fold;
 85:       VecCopy(blmP->Xold, tao->solution);
 86:       VecCopy(blmP->Gold, blmP->unprojected_gradient);

 88:       if (f != 0.0) {
 89:         delta = 2.0* PetscAbsScalar(f) / (gnorm*gnorm);
 90:       } else {
 91:         delta = 2.0/ (gnorm*gnorm);
 92:       }
 93:       MatLMVMSetDelta(blmP->M,delta);
 94:       MatLMVMReset(blmP->M);
 95:       MatLMVMUpdate(blmP->M, tao->solution, blmP->unprojected_gradient);
 96:       MatLMVMSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection);
 97:       VecScale(tao->stepdirection, -1.0);

 99:       /* This may be incorrect; linesearch has values fo stepmax and stepmin
100:          that should be reset. */
101:       TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
102:       TaoLineSearchApply(tao->linesearch,tao->solution,&f, blmP->unprojected_gradient, tao->stepdirection,  &stepsize, &ls_status);
103:       TaoAddLineSearchCounts(tao);

105:       if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
106:         tao->reason = TAO_DIVERGED_LS_FAILURE;
107:         break;
108:       }
109:     }

111:     /* Check for converged */
112:     VecBoundGradientProjection(blmP->unprojected_gradient, tao->solution, tao->XL, tao->XU, tao->gradient);
113:     VecNorm(tao->gradient, NORM_2, &gnorm);


116:     if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Not-a-Number");
117:     tao->niter++;
118:     TaoMonitor(tao, tao->niter, f, gnorm, 0.0, stepsize, &reason);
119:   }
120:   return(0);
121: }

125: static PetscErrorCode TaoSetup_BLMVM(Tao tao)
126: {
127:   TAO_BLMVM      *blmP = (TAO_BLMVM *)tao->data;
128:   PetscInt       n,N;

132:   /* Existence of tao->solution checked in TaoSetup() */
133:   VecDuplicate(tao->solution,&blmP->Xold);
134:   VecDuplicate(tao->solution,&blmP->Gold);
135:   VecDuplicate(tao->solution, &blmP->unprojected_gradient);

137:   if (!tao->stepdirection) {
138:     VecDuplicate(tao->solution, &tao->stepdirection);
139:   }
140:   if (!tao->gradient) {
141:     VecDuplicate(tao->solution,&tao->gradient);
142:   }
143:   if (!tao->XL) {
144:     VecDuplicate(tao->solution,&tao->XL);
145:     VecSet(tao->XL,PETSC_NINFINITY);
146:   }
147:   if (!tao->XU) {
148:     VecDuplicate(tao->solution,&tao->XU);
149:     VecSet(tao->XU,PETSC_INFINITY);
150:   }
151:   /* Create matrix for the limited memory approximation */
152:   VecGetLocalSize(tao->solution,&n);
153:   VecGetSize(tao->solution,&N);
154:   MatCreateLMVM(((PetscObject)tao)->comm,n,N,&blmP->M);
155:   MatLMVMAllocateVectors(blmP->M,tao->solution);
156:   return(0);
157: }

159: /* ---------------------------------------------------------- */
162: static PetscErrorCode TaoDestroy_BLMVM(Tao tao)
163: {
164:   TAO_BLMVM      *blmP = (TAO_BLMVM *)tao->data;

168:   if (tao->setupcalled) {
169:     MatDestroy(&blmP->M);
170:     VecDestroy(&blmP->unprojected_gradient);
171:     VecDestroy(&blmP->Xold);
172:     VecDestroy(&blmP->Gold);
173:   }
174:   PetscFree(tao->data);
175:   return(0);
176: }

178: /*------------------------------------------------------------*/
181: static PetscErrorCode TaoSetFromOptions_BLMVM(PetscOptions* PetscOptionsObject,Tao tao)
182: {

186:   PetscOptionsHead(PetscOptionsObject,"Limited-memory variable-metric method for bound constrained optimization");
187:   TaoLineSearchSetFromOptions(tao->linesearch);
188:   PetscOptionsTail();
189:   return(0);
190: }


193: /*------------------------------------------------------------*/
196: static int TaoView_BLMVM(Tao tao, PetscViewer viewer)
197: {
198:   TAO_BLMVM      *lmP = (TAO_BLMVM *)tao->data;
199:   PetscBool      isascii;

203:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
204:   if (isascii) {
205:     PetscViewerASCIIPushTab(viewer);
206:     PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", lmP->grad);
207:     PetscViewerASCIIPopTab(viewer);
208:   }
209:   return(0);
210: }

214: static PetscErrorCode TaoComputeDual_BLMVM(Tao tao, Vec DXL, Vec DXU)
215: {
216:   TAO_BLMVM      *blm = (TAO_BLMVM *) tao->data;

223:   if (!tao->gradient || !blm->unprojected_gradient) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Dual variables don't exist yet or no longer exist.\n");

225:   VecCopy(tao->gradient,DXL);
226:   VecAXPY(DXL,-1.0,blm->unprojected_gradient);
227:   VecSet(DXU,0.0);
228:   VecPointwiseMax(DXL,DXL,DXU);

230:   VecCopy(blm->unprojected_gradient,DXU);
231:   VecAXPY(DXU,-1.0,tao->gradient);
232:   VecAXPY(DXU,1.0,DXL);
233:   return(0);
234: }

236: /* ---------------------------------------------------------- */
237: /*MC
238:   TAOBLMVM - Bounded limited memory variable metric is a quasi-Newton method
239:          for nonlinear minimization with bound constraints. It is an extension
240:          of TAOLMVM

242:   Options Database Keys:
243: +     -tao_lmm_vectors - number of vectors to use for approximation
244: .     -tao_lmm_scale_type - "none","scalar","broyden"
245: .     -tao_lmm_limit_type - "none","average","relative","absolute"
246: .     -tao_lmm_rescale_type - "none","scalar","gl"
247: .     -tao_lmm_limit_mu - mu limiting factor
248: .     -tao_lmm_limit_nu - nu limiting factor
249: .     -tao_lmm_delta_min - minimum delta value
250: .     -tao_lmm_delta_max - maximum delta value
251: .     -tao_lmm_broyden_phi - phi factor for Broyden scaling
252: .     -tao_lmm_scalar_alpha - alpha factor for scalar scaling
253: .     -tao_lmm_rescale_alpha - alpha factor for rescaling diagonal
254: .     -tao_lmm_rescale_beta - beta factor for rescaling diagonal
255: .     -tao_lmm_scalar_history - amount of history for scalar scaling
256: .     -tao_lmm_rescale_history - amount of history for rescaling diagonal
257: -     -tao_lmm_eps - rejection tolerance

259:   Level: beginner
260: M*/
263: PETSC_EXTERN PetscErrorCode TaoCreate_BLMVM(Tao tao)
264: {
265:   TAO_BLMVM      *blmP;
266:   const char     *morethuente_type = TAOLINESEARCHMT;

270:   tao->ops->setup = TaoSetup_BLMVM;
271:   tao->ops->solve = TaoSolve_BLMVM;
272:   tao->ops->view = TaoView_BLMVM;
273:   tao->ops->setfromoptions = TaoSetFromOptions_BLMVM;
274:   tao->ops->destroy = TaoDestroy_BLMVM;
275:   tao->ops->computedual = TaoComputeDual_BLMVM;

277:   PetscNewLog(tao,&blmP);
278:   tao->data = (void*)blmP;

280:   /* Override default settings (unless already changed) */
281:   if (!tao->max_it_changed) tao->max_it = 2000;
282:   if (!tao->max_funcs_changed) tao->max_funcs = 4000;
283:   if (!tao->fatol_changed) tao->fatol = 1.0e-4;
284:   if (!tao->frtol_changed) tao->frtol = 1.0e-4;

286:   TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
287:   TaoLineSearchSetType(tao->linesearch, morethuente_type);
288:   TaoLineSearchUseTaoRoutines(tao->linesearch,tao);
289:   TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
290:   return(0);
291: }