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: PetscUseTypeMethod(tao, convergencetest, 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: PetscUseTypeMethod(tao, update, tao->niter, tao->user_update);
43: TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
44: }
45: /* Compute direction */
46: gnorm2 = gnorm * gnorm;
47: if (gnorm2 == 0.0) gnorm2 = PETSC_MACHINE_EPSILON;
48: if (f == 0.0) {
49: delta = 2.0 / gnorm2;
50: } else {
51: delta = 2.0 * PetscAbsScalar(f) / gnorm2;
52: }
53: MatLMVMSymBroydenSetDelta(blmP->M, delta);
54: MatLMVMUpdate(blmP->M, tao->solution, tao->gradient);
55: MatSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection);
56: VecBoundGradientProjection(tao->stepdirection, tao->solution, tao->XL, tao->XU, tao->gradient);
58: /* Check for success (descent direction) */
59: VecDot(blmP->unprojected_gradient, tao->gradient, &gdx);
60: if (gdx <= 0) {
61: /* Step is not descent or solve was not successful
62: Use steepest descent direction (scaled) */
63: ++blmP->grad;
65: MatLMVMReset(blmP->M, PETSC_FALSE);
66: MatLMVMUpdate(blmP->M, tao->solution, blmP->unprojected_gradient);
67: MatSolve(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: MatLMVMReset(blmP->M, PETSC_FALSE);
89: MatLMVMUpdate(blmP->M, tao->solution, blmP->unprojected_gradient);
90: MatSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection);
91: VecScale(tao->stepdirection, -1.0);
93: /* This may be incorrect; linesearch has values for stepmax and stepmin
94: that should be reset. */
95: TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0);
96: TaoLineSearchApply(tao->linesearch, tao->solution, &f, blmP->unprojected_gradient, tao->stepdirection, &stepsize, &ls_status);
97: TaoAddLineSearchCounts(tao);
99: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
100: tao->reason = TAO_DIVERGED_LS_FAILURE;
101: break;
102: }
103: }
105: /* Check for converged */
106: VecBoundGradientProjection(blmP->unprojected_gradient, tao->solution, tao->XL, tao->XU, tao->gradient);
107: TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm);
109: tao->niter++;
110: TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its);
111: TaoMonitor(tao, tao->niter, f, gnorm, 0.0, stepsize);
112: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
113: }
114: return 0;
115: }
117: static PetscErrorCode TaoSetup_BLMVM(Tao tao)
118: {
119: TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data;
121: /* Existence of tao->solution checked in TaoSetup() */
122: VecDuplicate(tao->solution, &blmP->Xold);
123: VecDuplicate(tao->solution, &blmP->Gold);
124: VecDuplicate(tao->solution, &blmP->unprojected_gradient);
125: if (!tao->stepdirection) VecDuplicate(tao->solution, &tao->stepdirection);
126: if (!tao->gradient) VecDuplicate(tao->solution, &tao->gradient);
127: /* Allocate matrix for the limited memory approximation */
128: MatLMVMAllocate(blmP->M, tao->solution, blmP->unprojected_gradient);
130: /* If the user has set a matrix to solve as the initial H0, set the options prefix here, and set up the KSP */
131: if (blmP->H0) MatLMVMSetJ0(blmP->M, blmP->H0);
132: return 0;
133: }
135: /* ---------------------------------------------------------- */
136: static PetscErrorCode TaoDestroy_BLMVM(Tao tao)
137: {
138: TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data;
140: if (tao->setupcalled) {
141: VecDestroy(&blmP->unprojected_gradient);
142: VecDestroy(&blmP->Xold);
143: VecDestroy(&blmP->Gold);
144: }
145: MatDestroy(&blmP->M);
146: if (blmP->H0) PetscObjectDereference((PetscObject)blmP->H0);
147: PetscFree(tao->data);
148: return 0;
149: }
151: /*------------------------------------------------------------*/
152: static PetscErrorCode TaoSetFromOptions_BLMVM(Tao tao, PetscOptionItems *PetscOptionsObject)
153: {
154: TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data;
155: PetscBool is_spd, is_set;
157: PetscOptionsHeadBegin(PetscOptionsObject, "Limited-memory variable-metric method for bound constrained optimization");
158: PetscOptionsBool("-tao_blmvm_recycle", "enable recycling of the BFGS matrix between subsequent TaoSolve() calls", "", blmP->recycle, &blmP->recycle, NULL);
159: PetscOptionsHeadEnd();
160: MatSetOptionsPrefix(blmP->M, ((PetscObject)tao)->prefix);
161: MatAppendOptionsPrefix(blmP->M, "tao_blmvm_");
162: MatSetFromOptions(blmP->M);
163: MatIsSPDKnown(blmP->M, &is_set, &is_spd);
165: return 0;
166: }
168: /*------------------------------------------------------------*/
169: static PetscErrorCode TaoView_BLMVM(Tao tao, PetscViewer viewer)
170: {
171: TAO_BLMVM *lmP = (TAO_BLMVM *)tao->data;
172: PetscBool isascii;
174: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
175: if (isascii) {
176: PetscViewerASCIIPrintf(viewer, "Gradient steps: %" PetscInt_FMT "\n", lmP->grad);
177: PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO);
178: MatView(lmP->M, viewer);
179: PetscViewerPopFormat(viewer);
180: }
181: return 0;
182: }
184: static PetscErrorCode TaoComputeDual_BLMVM(Tao tao, Vec DXL, Vec DXU)
185: {
186: TAO_BLMVM *blm = (TAO_BLMVM *)tao->data;
193: VecCopy(tao->gradient, DXL);
194: VecAXPY(DXL, -1.0, blm->unprojected_gradient);
195: VecSet(DXU, 0.0);
196: VecPointwiseMax(DXL, DXL, DXU);
198: VecCopy(blm->unprojected_gradient, DXU);
199: VecAXPY(DXU, -1.0, tao->gradient);
200: VecAXPY(DXU, 1.0, DXL);
201: return 0;
202: }
204: /* ---------------------------------------------------------- */
205: /*MC
206: TAOBLMVM - Bounded limited memory variable metric is a quasi-Newton method
207: for nonlinear minimization with bound constraints. It is an extension
208: of TAOLMVM
210: Options Database Keys:
211: . -tao_lmm_recycle - enable recycling of LMVM information between subsequent TaoSolve calls
213: Level: beginner
214: M*/
215: PETSC_EXTERN PetscErrorCode TaoCreate_BLMVM(Tao tao)
216: {
217: TAO_BLMVM *blmP;
218: const char *morethuente_type = TAOLINESEARCHMT;
220: tao->ops->setup = TaoSetup_BLMVM;
221: tao->ops->solve = TaoSolve_BLMVM;
222: tao->ops->view = TaoView_BLMVM;
223: tao->ops->setfromoptions = TaoSetFromOptions_BLMVM;
224: tao->ops->destroy = TaoDestroy_BLMVM;
225: tao->ops->computedual = TaoComputeDual_BLMVM;
227: PetscNew(&blmP);
228: blmP->H0 = NULL;
229: blmP->recycle = PETSC_FALSE;
230: tao->data = (void *)blmP;
232: /* Override default settings (unless already changed) */
233: if (!tao->max_it_changed) tao->max_it = 2000;
234: if (!tao->max_funcs_changed) tao->max_funcs = 4000;
236: TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
237: PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
238: TaoLineSearchSetType(tao->linesearch, morethuente_type);
239: TaoLineSearchUseTaoRoutines(tao->linesearch, tao);
241: KSPInitializePackage();
242: MatCreate(((PetscObject)tao)->comm, &blmP->M);
243: MatSetType(blmP->M, MATLMVMBFGS);
244: PetscObjectIncrementTabLevel((PetscObject)blmP->M, (PetscObject)tao, 1);
245: return 0;
246: }
248: /*@
249: TaoLMVMRecycle - Enable/disable recycling of the QN history between subsequent TaoSolve calls.
251: Input Parameters:
252: + tao - the Tao solver context
253: - flg - Boolean flag for recycling (PETSC_TRUE or PETSC_FALSE)
255: Level: intermediate
256: @*/
257: PetscErrorCode TaoLMVMRecycle(Tao tao, PetscBool flg)
258: {
259: TAO_LMVM *lmP;
260: TAO_BLMVM *blmP;
261: PetscBool is_lmvm, is_blmvm;
263: PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm);
264: PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm);
265: if (is_lmvm) {
266: lmP = (TAO_LMVM *)tao->data;
267: lmP->recycle = flg;
268: } else if (is_blmvm) {
269: blmP = (TAO_BLMVM *)tao->data;
270: blmP->recycle = flg;
271: }
272: return 0;
273: }
275: /*@
276: TaoLMVMSetH0 - Set the initial Hessian for the QN approximation
278: Input Parameters:
279: + tao - the Tao solver context
280: - H0 - Mat object for the initial Hessian
282: Level: advanced
284: .seealso: `TaoLMVMGetH0()`, `TaoLMVMGetH0KSP()`
285: @*/
286: PetscErrorCode TaoLMVMSetH0(Tao tao, Mat H0)
287: {
288: TAO_LMVM *lmP;
289: TAO_BLMVM *blmP;
290: PetscBool is_lmvm, is_blmvm;
292: PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm);
293: PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm);
294: if (is_lmvm) {
295: lmP = (TAO_LMVM *)tao->data;
296: PetscObjectReference((PetscObject)H0);
297: lmP->H0 = H0;
298: } else if (is_blmvm) {
299: blmP = (TAO_BLMVM *)tao->data;
300: PetscObjectReference((PetscObject)H0);
301: blmP->H0 = H0;
302: }
303: return 0;
304: }
306: /*@
307: TaoLMVMGetH0 - Get the matrix object for the QN initial Hessian
309: Input Parameters:
310: . tao - the Tao solver context
312: Output Parameters:
313: . H0 - Mat object for the initial Hessian
315: Level: advanced
317: .seealso: `TaoLMVMSetH0()`, `TaoLMVMGetH0KSP()`
318: @*/
319: PetscErrorCode TaoLMVMGetH0(Tao tao, Mat *H0)
320: {
321: TAO_LMVM *lmP;
322: TAO_BLMVM *blmP;
323: PetscBool is_lmvm, is_blmvm;
324: Mat M;
326: PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm);
327: PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm);
328: if (is_lmvm) {
329: lmP = (TAO_LMVM *)tao->data;
330: M = lmP->M;
331: } else if (is_blmvm) {
332: blmP = (TAO_BLMVM *)tao->data;
333: M = blmP->M;
334: } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONG, "This routine applies to TAO_LMVM and TAO_BLMVM.");
335: MatLMVMGetJ0(M, H0);
336: return 0;
337: }
339: /*@
340: TaoLMVMGetH0KSP - Get the iterative solver for applying the inverse of the QN initial Hessian
342: Input Parameters:
343: . tao - the Tao solver context
345: Output Parameters:
346: . ksp - KSP solver context for the initial Hessian
348: Level: advanced
350: .seealso: `TaoLMVMGetH0()`, `TaoLMVMGetH0KSP()`
351: @*/
352: PetscErrorCode TaoLMVMGetH0KSP(Tao tao, KSP *ksp)
353: {
354: TAO_LMVM *lmP;
355: TAO_BLMVM *blmP;
356: PetscBool is_lmvm, is_blmvm;
357: Mat M;
359: PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm);
360: PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm);
361: if (is_lmvm) {
362: lmP = (TAO_LMVM *)tao->data;
363: M = lmP->M;
364: } else if (is_blmvm) {
365: blmP = (TAO_BLMVM *)tao->data;
366: M = blmP->M;
367: } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONG, "This routine applies to TAO_LMVM and TAO_BLMVM.");
368: MatLMVMGetJ0KSP(M, ksp);
369: return 0;
370: }