Actual source code: blmvm.c

petsc-3.7.7 2017-09-25
Report Typos and Errors
  1: #include <petsctaolinesearch.h>
  2: #include <../src/tao/matrix/lmvmmat.h>
  3: #include <../src/tao/unconstrained/impls/lmvm/lmvm.h>
  4: #include <../src/tao/bound/impls/blmvm/blmvm.h>

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

127: static PetscErrorCode TaoSetup_BLMVM(Tao tao)
128: {
129:   TAO_BLMVM      *blmP = (TAO_BLMVM *)tao->data;
130:   PetscInt       n,N;
132:   KSP            H0ksp;

135:   /* Existence of tao->solution checked in TaoSetup() */
136:   VecDuplicate(tao->solution,&blmP->Xold);
137:   VecDuplicate(tao->solution,&blmP->Gold);
138:   VecDuplicate(tao->solution, &blmP->unprojected_gradient);

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

160:   /* If the user has set a matrix to solve as the initial H0, set the options prefix here, and set up the KSP */
161:   if (blmP->H0) {
162:     const char *prefix;
163:     PC H0pc;

165:     MatLMVMSetH0(blmP->M, blmP->H0);
166:     MatLMVMGetH0KSP(blmP->M, &H0ksp);

168:     TaoGetOptionsPrefix(tao, &prefix);
169:     KSPSetOptionsPrefix(H0ksp, prefix);
170:     PetscObjectAppendOptionsPrefix((PetscObject)H0ksp, "tao_h0_");
171:     KSPGetPC(H0ksp, &H0pc);
172:     PetscObjectAppendOptionsPrefix((PetscObject)H0pc,  "tao_h0_");

174:     KSPSetFromOptions(H0ksp);
175:     KSPSetUp(H0ksp);
176:   }

178:   return(0);
179: }

181: /* ---------------------------------------------------------- */
184: static PetscErrorCode TaoDestroy_BLMVM(Tao tao)
185: {
186:   TAO_BLMVM      *blmP = (TAO_BLMVM *)tao->data;

190:   if (tao->setupcalled) {
191:     MatDestroy(&blmP->M);
192:     VecDestroy(&blmP->unprojected_gradient);
193:     VecDestroy(&blmP->Xold);
194:     VecDestroy(&blmP->Gold);
195:   }

197:   if (blmP->H0) {
198:     PetscObjectDereference((PetscObject)blmP->H0);
199:   }

201:   PetscFree(tao->data);
202:   return(0);
203: }

205: /*------------------------------------------------------------*/
208: static PetscErrorCode TaoSetFromOptions_BLMVM(PetscOptionItems* PetscOptionsObject,Tao tao)
209: {

213:   PetscOptionsHead(PetscOptionsObject,"Limited-memory variable-metric method for bound constrained optimization");
214:   TaoLineSearchSetFromOptions(tao->linesearch);
215:   PetscOptionsTail();
216:   return(0);
217: }


220: /*------------------------------------------------------------*/
223: static int TaoView_BLMVM(Tao tao, PetscViewer viewer)
224: {
225:   TAO_BLMVM      *lmP = (TAO_BLMVM *)tao->data;
226:   PetscBool      isascii;

230:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
231:   if (isascii) {
232:     PetscViewerASCIIPushTab(viewer);
233:     PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", lmP->grad);
234:     PetscViewerASCIIPopTab(viewer);
235:   }
236:   return(0);
237: }

241: static PetscErrorCode TaoComputeDual_BLMVM(Tao tao, Vec DXL, Vec DXU)
242: {
243:   TAO_BLMVM      *blm = (TAO_BLMVM *) tao->data;

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

252:   VecCopy(tao->gradient,DXL);
253:   VecAXPY(DXL,-1.0,blm->unprojected_gradient);
254:   VecSet(DXU,0.0);
255:   VecPointwiseMax(DXL,DXL,DXU);

257:   VecCopy(blm->unprojected_gradient,DXU);
258:   VecAXPY(DXU,-1.0,tao->gradient);
259:   VecAXPY(DXU,1.0,DXL);
260:   return(0);
261: }

263: /* ---------------------------------------------------------- */
264: /*MC
265:   TAOBLMVM - Bounded limited memory variable metric is a quasi-Newton method
266:          for nonlinear minimization with bound constraints. It is an extension
267:          of TAOLMVM

269:   Options Database Keys:
270: +     -tao_lmm_vectors - number of vectors to use for approximation
271: .     -tao_lmm_scale_type - "none","scalar","broyden"
272: .     -tao_lmm_limit_type - "none","average","relative","absolute"
273: .     -tao_lmm_rescale_type - "none","scalar","gl"
274: .     -tao_lmm_limit_mu - mu limiting factor
275: .     -tao_lmm_limit_nu - nu limiting factor
276: .     -tao_lmm_delta_min - minimum delta value
277: .     -tao_lmm_delta_max - maximum delta value
278: .     -tao_lmm_broyden_phi - phi factor for Broyden scaling
279: .     -tao_lmm_scalar_alpha - alpha factor for scalar scaling
280: .     -tao_lmm_rescale_alpha - alpha factor for rescaling diagonal
281: .     -tao_lmm_rescale_beta - beta factor for rescaling diagonal
282: .     -tao_lmm_scalar_history - amount of history for scalar scaling
283: .     -tao_lmm_rescale_history - amount of history for rescaling diagonal
284: -     -tao_lmm_eps - rejection tolerance

286:   Level: beginner
287: M*/
290: PETSC_EXTERN PetscErrorCode TaoCreate_BLMVM(Tao tao)
291: {
292:   TAO_BLMVM      *blmP;
293:   const char     *morethuente_type = TAOLINESEARCHMT;

297:   tao->ops->setup = TaoSetup_BLMVM;
298:   tao->ops->solve = TaoSolve_BLMVM;
299:   tao->ops->view = TaoView_BLMVM;
300:   tao->ops->setfromoptions = TaoSetFromOptions_BLMVM;
301:   tao->ops->destroy = TaoDestroy_BLMVM;
302:   tao->ops->computedual = TaoComputeDual_BLMVM;

304:   PetscNewLog(tao,&blmP);
305:   blmP->H0 = NULL;
306:   tao->data = (void*)blmP;

308:   /* Override default settings (unless already changed) */
309:   if (!tao->max_it_changed) tao->max_it = 2000;
310:   if (!tao->max_funcs_changed) tao->max_funcs = 4000;

312:   TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
313:   TaoLineSearchSetType(tao->linesearch, morethuente_type);
314:   TaoLineSearchUseTaoRoutines(tao->linesearch,tao);
315:   TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
316:   return(0);
317: }

321: PETSC_EXTERN PetscErrorCode TaoLMVMSetH0(Tao tao, Mat H0)
322: {
323:   TAO_LMVM       *lmP;
324:   TAO_BLMVM      *blmP;
325:   const TaoType  type;
326:   PetscBool is_lmvm, is_blmvm;


330:   TaoGetType(tao, &type);
331:   PetscStrcmp(type, TAOLMVM,  &is_lmvm);
332:   PetscStrcmp(type, TAOBLMVM, &is_blmvm);

334:   if (is_lmvm) {
335:     lmP = (TAO_LMVM *)tao->data;
336:     PetscObjectReference((PetscObject)H0);
337:     lmP->H0 = H0;
338:   } else if (is_blmvm) {
339:     blmP = (TAO_BLMVM *)tao->data;
340:     PetscObjectReference((PetscObject)H0);
341:     blmP->H0 = H0;
342:   } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONGSTATE, "This routine applies to TAO_LMVM and TAO_BLMVM.");

344:   return(0);
345: }

349: PETSC_EXTERN PetscErrorCode TaoLMVMGetH0(Tao tao, Mat *H0)
350: {
351:   TAO_LMVM       *lmP;
352:   TAO_BLMVM      *blmP;
353:   const TaoType  type;
354:   PetscBool      is_lmvm, is_blmvm;
355:   Mat            M;


359:   TaoGetType(tao, &type);
360:   PetscStrcmp(type, TAOLMVM,  &is_lmvm);
361:   PetscStrcmp(type, TAOBLMVM, &is_blmvm);

363:   if (is_lmvm) {
364:     lmP = (TAO_LMVM *)tao->data;
365:     M = lmP->M;
366:   } else if (is_blmvm) {
367:     blmP = (TAO_BLMVM *)tao->data;
368:     M = blmP->M;
369:   } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONGSTATE, "This routine applies to TAO_LMVM and TAO_BLMVM.");

371:   MatLMVMGetH0(M, H0);
372:   return(0);
373: }

377: PETSC_EXTERN PetscErrorCode TaoLMVMGetH0KSP(Tao tao, KSP *ksp)
378: {
379:   TAO_LMVM       *lmP;
380:   TAO_BLMVM      *blmP;
381:   const TaoType  type;
382:   PetscBool      is_lmvm, is_blmvm;
383:   Mat            M;

386:   TaoGetType(tao, &type);
387:   PetscStrcmp(type, TAOLMVM,  &is_lmvm);
388:   PetscStrcmp(type, TAOBLMVM, &is_blmvm);

390:   if (is_lmvm) {
391:     lmP = (TAO_LMVM *)tao->data;
392:     M = lmP->M;
393:   } else if (is_blmvm) {
394:     blmP = (TAO_BLMVM *)tao->data;
395:     M = blmP->M;
396:   } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONGSTATE, "This routine applies to TAO_LMVM and TAO_BLMVM.");

398:   MatLMVMGetH0KSP(M, ksp);
399:   return(0);
400: }