Actual source code: blmvm.c

  1: #include <petsctaolinesearch.h>
  2: #include <../src/tao/unconstrained/impls/lmvm/lmvm.h>
  3: #include <../src/tao/bound/impls/blmvm/blmvm.h>

  5: /*------------------------------------------------------------*/
  6: static PetscErrorCode TaoSolve_BLMVM(Tao tao)
  7: {
  8:   TAO_BLMVM                    *blmP = (TAO_BLMVM *)tao->data;
  9:   TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;
 10:   PetscReal                    f, fold, gdx, gnorm, gnorm2;
 11:   PetscReal                    stepsize = 1.0,delta;

 13:   /*  Project initial point onto bounds */
 14:   TaoComputeVariableBounds(tao);
 15:   VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);
 16:   TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);

 18:   /* Check convergence criteria */
 19:   TaoComputeObjectiveAndGradient(tao, tao->solution,&f,blmP->unprojected_gradient);
 20:   VecBoundGradientProjection(blmP->unprojected_gradient,tao->solution, tao->XL,tao->XU,tao->gradient);

 22:   TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);

 25:   tao->reason = TAO_CONTINUE_ITERATING;
 26:   TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
 27:   TaoMonitor(tao,tao->niter,f,gnorm,0.0,stepsize);
 28:   (*tao->ops->convergencetest)(tao,tao->cnvP);
 29:   if (tao->reason != TAO_CONTINUE_ITERATING) return 0;

 31:   /* Set counter for gradient/reset steps */
 32:   if (!blmP->recycle) {
 33:     blmP->grad = 0;
 34:     blmP->reset = 0;
 35:     MatLMVMReset(blmP->M, PETSC_FALSE);
 36:   }

 38:   /* Have not converged; continue with Newton method */
 39:   while (tao->reason == TAO_CONTINUE_ITERATING) {
 40:     /* Call general purpose update function */
 41:     if (tao->ops->update) {
 42:       (*tao->ops->update)(tao, tao->niter, tao->user_update);
 43:     }
 44:     /* Compute direction */
 45:     gnorm2 = gnorm*gnorm;
 46:     if (gnorm2 == 0.0) gnorm2 = PETSC_MACHINE_EPSILON;
 47:     if (f == 0.0) {
 48:       delta = 2.0 / gnorm2;
 49:     } else {
 50:       delta = 2.0 * PetscAbsScalar(f) / gnorm2;
 51:     }
 52:     MatLMVMSymBroydenSetDelta(blmP->M, delta);
 53:     MatLMVMUpdate(blmP->M, tao->solution, tao->gradient);
 54:     MatSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection);
 55:     VecBoundGradientProjection(tao->stepdirection,tao->solution,tao->XL,tao->XU,tao->gradient);

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

 64:       MatLMVMReset(blmP->M, PETSC_FALSE);
 65:       MatLMVMUpdate(blmP->M, tao->solution, blmP->unprojected_gradient);
 66:       MatSolve(blmP->M,blmP->unprojected_gradient, tao->stepdirection);
 67:     }
 68:     VecScale(tao->stepdirection,-1.0);

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

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

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

 87:       MatLMVMReset(blmP->M, PETSC_FALSE);
 88:       MatLMVMUpdate(blmP->M, tao->solution, blmP->unprojected_gradient);
 89:       MatSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection);
 90:       VecScale(tao->stepdirection, -1.0);

 92:       /* This may be incorrect; linesearch has values for stepmax and stepmin
 93:          that should be reset. */
 94:       TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
 95:       TaoLineSearchApply(tao->linesearch,tao->solution,&f, blmP->unprojected_gradient, tao->stepdirection,  &stepsize, &ls_status);
 96:       TaoAddLineSearchCounts(tao);

 98:       if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
 99:         tao->reason = TAO_DIVERGED_LS_FAILURE;
100:         break;
101:       }
102:     }

104:     /* Check for converged */
105:     VecBoundGradientProjection(blmP->unprojected_gradient, tao->solution, tao->XL, tao->XU, tao->gradient);
106:     TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm);
108:     tao->niter++;
109:     TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
110:     TaoMonitor(tao,tao->niter,f,gnorm,0.0,stepsize);
111:     (*tao->ops->convergencetest)(tao,tao->cnvP);
112:   }
113:   return 0;
114: }

116: static PetscErrorCode TaoSetup_BLMVM(Tao tao)
117: {
118:   TAO_BLMVM      *blmP = (TAO_BLMVM *)tao->data;

120:   /* Existence of tao->solution checked in TaoSetup() */
121:   VecDuplicate(tao->solution,&blmP->Xold);
122:   VecDuplicate(tao->solution,&blmP->Gold);
123:   VecDuplicate(tao->solution, &blmP->unprojected_gradient);

125:   if (!tao->stepdirection) {
126:     VecDuplicate(tao->solution, &tao->stepdirection);
127:   }
128:   if (!tao->gradient) {
129:     VecDuplicate(tao->solution,&tao->gradient);
130:   }
131:   if (!tao->XL) {
132:     VecDuplicate(tao->solution,&tao->XL);
133:     VecSet(tao->XL,PETSC_NINFINITY);
134:   }
135:   if (!tao->XU) {
136:     VecDuplicate(tao->solution,&tao->XU);
137:     VecSet(tao->XU,PETSC_INFINITY);
138:   }
139:   /* Allocate matrix for the limited memory approximation */
140:   MatLMVMAllocate(blmP->M,tao->solution,blmP->unprojected_gradient);

142:   /* If the user has set a matrix to solve as the initial H0, set the options prefix here, and set up the KSP */
143:   if (blmP->H0) {
144:     MatLMVMSetJ0(blmP->M, blmP->H0);
145:   }
146:   return 0;
147: }

149: /* ---------------------------------------------------------- */
150: static PetscErrorCode TaoDestroy_BLMVM(Tao tao)
151: {
152:   TAO_BLMVM      *blmP = (TAO_BLMVM *)tao->data;

154:   if (tao->setupcalled) {
155:     VecDestroy(&blmP->unprojected_gradient);
156:     VecDestroy(&blmP->Xold);
157:     VecDestroy(&blmP->Gold);
158:   }
159:   MatDestroy(&blmP->M);
160:   if (blmP->H0) {
161:     PetscObjectDereference((PetscObject)blmP->H0);
162:   }
163:   PetscFree(tao->data);
164:   return 0;
165: }

167: /*------------------------------------------------------------*/
168: static PetscErrorCode TaoSetFromOptions_BLMVM(PetscOptionItems* PetscOptionsObject,Tao tao)
169: {
170:   TAO_BLMVM      *blmP = (TAO_BLMVM *)tao->data;
171:   PetscBool      is_spd;

173:   PetscOptionsHead(PetscOptionsObject,"Limited-memory variable-metric method for bound constrained optimization");
174:   PetscOptionsBool("-tao_blmvm_recycle","enable recycling of the BFGS matrix between subsequent TaoSolve() calls","",blmP->recycle,&blmP->recycle,NULL);
175:   PetscOptionsTail();
176:   MatSetOptionsPrefix(blmP->M, ((PetscObject)tao)->prefix);
177:   MatAppendOptionsPrefix(blmP->M, "tao_blmvm_");
178:   MatSetFromOptions(blmP->M);
179:   MatGetOption(blmP->M, MAT_SPD, &is_spd);
181:   return 0;
182: }

184: /*------------------------------------------------------------*/
185: static PetscErrorCode TaoView_BLMVM(Tao tao, PetscViewer viewer)
186: {
187:   TAO_BLMVM      *lmP = (TAO_BLMVM *)tao->data;
188:   PetscBool      isascii;

190:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
191:   if (isascii) {
192:     PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", lmP->grad);
193:     PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO);
194:     MatView(lmP->M, viewer);
195:     PetscViewerPopFormat(viewer);
196:   }
197:   return 0;
198: }

200: static PetscErrorCode TaoComputeDual_BLMVM(Tao tao, Vec DXL, Vec DXU)
201: {
202:   TAO_BLMVM      *blm = (TAO_BLMVM *) tao->data;


209:   VecCopy(tao->gradient,DXL);
210:   VecAXPY(DXL,-1.0,blm->unprojected_gradient);
211:   VecSet(DXU,0.0);
212:   VecPointwiseMax(DXL,DXL,DXU);

214:   VecCopy(blm->unprojected_gradient,DXU);
215:   VecAXPY(DXU,-1.0,tao->gradient);
216:   VecAXPY(DXU,1.0,DXL);
217:   return 0;
218: }

220: /* ---------------------------------------------------------- */
221: /*MC
222:   TAOBLMVM - Bounded limited memory variable metric is a quasi-Newton method
223:          for nonlinear minimization with bound constraints. It is an extension
224:          of TAOLMVM

226:   Options Database Keys:
227: .     -tao_lmm_recycle - enable recycling of LMVM information between subsequent TaoSolve calls

229:   Level: beginner
230: M*/
231: PETSC_EXTERN PetscErrorCode TaoCreate_BLMVM(Tao tao)
232: {
233:   TAO_BLMVM      *blmP;
234:   const char     *morethuente_type = TAOLINESEARCHMT;

236:   tao->ops->setup = TaoSetup_BLMVM;
237:   tao->ops->solve = TaoSolve_BLMVM;
238:   tao->ops->view = TaoView_BLMVM;
239:   tao->ops->setfromoptions = TaoSetFromOptions_BLMVM;
240:   tao->ops->destroy = TaoDestroy_BLMVM;
241:   tao->ops->computedual = TaoComputeDual_BLMVM;

243:   PetscNewLog(tao,&blmP);
244:   blmP->H0 = NULL;
245:   blmP->recycle = PETSC_FALSE;
246:   tao->data = (void*)blmP;

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

252:   TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
253:   PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
254:   TaoLineSearchSetType(tao->linesearch, morethuente_type);
255:   TaoLineSearchUseTaoRoutines(tao->linesearch,tao);

257:   KSPInitializePackage();
258:   MatCreate(((PetscObject)tao)->comm, &blmP->M);
259:   MatSetType(blmP->M, MATLMVMBFGS);
260:   PetscObjectIncrementTabLevel((PetscObject)blmP->M, (PetscObject)tao, 1);
261:   return 0;
262: }

264: /*@
265:   TaoLMVMRecycle - Enable/disable recycling of the QN history between subsequent TaoSolve calls.

267:   Input Parameters:
268: +  tao  - the Tao solver context
269: -  flg - Boolean flag for recycling (PETSC_TRUE or PETSC_FALSE)

271:   Level: intermediate
272: @*/
273: PetscErrorCode TaoLMVMRecycle(Tao tao, PetscBool flg)
274: {
275:   TAO_LMVM       *lmP;
276:   TAO_BLMVM      *blmP;
277:   PetscBool      is_lmvm, is_blmvm;

279:   PetscObjectTypeCompare((PetscObject)tao,TAOLMVM,&is_lmvm);
280:   PetscObjectTypeCompare((PetscObject)tao,TAOBLMVM,&is_blmvm);
281:   if (is_lmvm) {
282:     lmP = (TAO_LMVM *)tao->data;
283:     lmP->recycle = flg;
284:   } else if (is_blmvm) {
285:     blmP = (TAO_BLMVM *)tao->data;
286:     blmP->recycle = flg;
287:   }
288:   return 0;
289: }

291: /*@
292:   TaoLMVMSetH0 - Set the initial Hessian for the QN approximation

294:   Input Parameters:
295: +  tao  - the Tao solver context
296: -  H0 - Mat object for the initial Hessian

298:   Level: advanced

300: .seealso: TaoLMVMGetH0(), TaoLMVMGetH0KSP()
301: @*/
302: PetscErrorCode TaoLMVMSetH0(Tao tao, Mat H0)
303: {
304:   TAO_LMVM       *lmP;
305:   TAO_BLMVM      *blmP;
306:   PetscBool      is_lmvm, is_blmvm;

308:   PetscObjectTypeCompare((PetscObject)tao,TAOLMVM,&is_lmvm);
309:   PetscObjectTypeCompare((PetscObject)tao,TAOBLMVM,&is_blmvm);
310:   if (is_lmvm) {
311:     lmP = (TAO_LMVM *)tao->data;
312:     PetscObjectReference((PetscObject)H0);
313:     lmP->H0 = H0;
314:   } else if (is_blmvm) {
315:     blmP = (TAO_BLMVM *)tao->data;
316:     PetscObjectReference((PetscObject)H0);
317:     blmP->H0 = H0;
318:   }
319:   return 0;
320: }

322: /*@
323:   TaoLMVMGetH0 - Get the matrix object for the QN initial Hessian

325:   Input Parameters:
326: .  tao  - the Tao solver context

328:   Output Parameters:
329: .  H0 - Mat object for the initial Hessian

331:   Level: advanced

333: .seealso: TaoLMVMSetH0(), TaoLMVMGetH0KSP()
334: @*/
335: PetscErrorCode TaoLMVMGetH0(Tao tao, Mat *H0)
336: {
337:   TAO_LMVM       *lmP;
338:   TAO_BLMVM      *blmP;
339:   PetscBool      is_lmvm, is_blmvm;
340:   Mat            M;

342:   PetscObjectTypeCompare((PetscObject)tao,TAOLMVM,&is_lmvm);
343:   PetscObjectTypeCompare((PetscObject)tao,TAOBLMVM,&is_blmvm);
344:   if (is_lmvm) {
345:     lmP = (TAO_LMVM *)tao->data;
346:     M = lmP->M;
347:   } else if (is_blmvm) {
348:     blmP = (TAO_BLMVM *)tao->data;
349:     M = blmP->M;
350:   } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONG, "This routine applies to TAO_LMVM and TAO_BLMVM.");
351:   MatLMVMGetJ0(M, H0);
352:   return 0;
353: }

355: /*@
356:   TaoLMVMGetH0KSP - Get the iterative solver for applying the inverse of the QN initial Hessian

358:   Input Parameters:
359: .  tao  - the Tao solver context

361:   Output Parameters:
362: .  ksp - KSP solver context for the initial Hessian

364:   Level: advanced

366: .seealso: TaoLMVMGetH0(), TaoLMVMGetH0KSP()
367: @*/
368: PetscErrorCode TaoLMVMGetH0KSP(Tao tao, KSP *ksp)
369: {
370:   TAO_LMVM       *lmP;
371:   TAO_BLMVM      *blmP;
372:   PetscBool      is_lmvm, is_blmvm;
373:   Mat            M;

375:   PetscObjectTypeCompare((PetscObject)tao,TAOLMVM,&is_lmvm);
376:   PetscObjectTypeCompare((PetscObject)tao,TAOBLMVM,&is_blmvm);
377:   if (is_lmvm) {
378:     lmP = (TAO_LMVM *)tao->data;
379:     M = lmP->M;
380:   } else if (is_blmvm) {
381:     blmP = (TAO_BLMVM *)tao->data;
382:     M = blmP->M;
383:   } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONG, "This routine applies to TAO_LMVM and TAO_BLMVM.");
384:   MatLMVMGetJ0KSP(M, ksp);
385:   return 0;
386: }