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