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: static PetscErrorCode TaoSolve_BLMVM(Tao tao)
6: {
7: TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data;
8: TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;
9: PetscReal f, fold, gdx, gnorm, gnorm2;
10: PetscReal stepsize = 1.0, delta;
12: PetscFunctionBegin;
13: /* Project initial point onto bounds */
14: PetscCall(TaoComputeVariableBounds(tao));
15: PetscCall(VecMedian(tao->XL, tao->solution, tao->XU, tao->solution));
16: PetscCall(TaoLineSearchSetVariableBounds(tao->linesearch, tao->XL, tao->XU));
18: /* Check convergence criteria */
19: PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, blmP->unprojected_gradient));
20: PetscCall(VecBoundGradientProjection(blmP->unprojected_gradient, tao->solution, tao->XL, tao->XU, tao->gradient));
22: PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm));
23: PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated infinity or NaN");
25: tao->reason = TAO_CONTINUE_ITERATING;
26: PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
27: PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, stepsize));
28: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
29: if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(PETSC_SUCCESS);
31: /* Set counter for gradient/reset steps */
32: if (!blmP->recycle) {
33: blmP->grad = 0;
34: blmP->reset = 0;
35: PetscCall(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: PetscCall(TaoComputeObjective(tao, tao->solution, &f));
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: PetscCall(MatLMVMSymBroydenSetDelta(blmP->M, delta));
54: PetscCall(MatLMVMUpdate(blmP->M, tao->solution, tao->gradient));
55: PetscCall(MatSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection));
56: PetscCall(VecBoundGradientProjection(tao->stepdirection, tao->solution, tao->XL, tao->XU, tao->gradient));
58: /* Check for success (descent direction) */
59: PetscCall(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: PetscCall(MatLMVMReset(blmP->M, PETSC_FALSE));
66: PetscCall(MatLMVMUpdate(blmP->M, tao->solution, blmP->unprojected_gradient));
67: PetscCall(MatSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection));
68: }
69: PetscCall(VecScale(tao->stepdirection, -1.0));
71: /* Perform the linesearch */
72: fold = f;
73: PetscCall(VecCopy(tao->solution, blmP->Xold));
74: PetscCall(VecCopy(blmP->unprojected_gradient, blmP->Gold));
75: PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0));
76: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, blmP->unprojected_gradient, tao->stepdirection, &stepsize, &ls_status));
77: PetscCall(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: PetscCall(VecCopy(blmP->Xold, tao->solution));
86: PetscCall(VecCopy(blmP->Gold, blmP->unprojected_gradient));
88: PetscCall(MatLMVMReset(blmP->M, PETSC_FALSE));
89: PetscCall(MatLMVMUpdate(blmP->M, tao->solution, blmP->unprojected_gradient));
90: PetscCall(MatSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection));
91: PetscCall(VecScale(tao->stepdirection, -1.0));
93: /* This may be incorrect; linesearch has values for stepmax and stepmin
94: that should be reset. */
95: PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0));
96: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, blmP->unprojected_gradient, tao->stepdirection, &stepsize, &ls_status));
97: PetscCall(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: PetscCall(VecBoundGradientProjection(blmP->unprojected_gradient, tao->solution, tao->XL, tao->XU, tao->gradient));
107: PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm));
108: PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Not-a-Number");
109: tao->niter++;
110: PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
111: PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, stepsize));
112: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
113: }
114: PetscFunctionReturn(PETSC_SUCCESS);
115: }
117: static PetscErrorCode TaoSetup_BLMVM(Tao tao)
118: {
119: TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data;
121: PetscFunctionBegin;
122: /* Existence of tao->solution checked in TaoSetup() */
123: PetscCall(VecDuplicate(tao->solution, &blmP->Xold));
124: PetscCall(VecDuplicate(tao->solution, &blmP->Gold));
125: PetscCall(VecDuplicate(tao->solution, &blmP->unprojected_gradient));
126: if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
127: if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient));
128: /* Allocate matrix for the limited memory approximation */
129: PetscCall(MatLMVMAllocate(blmP->M, tao->solution, blmP->unprojected_gradient));
131: /* If the user has set a matrix to solve as the initial H0, set the options prefix here, and set up the KSP */
132: if (blmP->H0) PetscCall(MatLMVMSetJ0(blmP->M, blmP->H0));
133: PetscFunctionReturn(PETSC_SUCCESS);
134: }
136: static PetscErrorCode TaoDestroy_BLMVM(Tao tao)
137: {
138: TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data;
140: PetscFunctionBegin;
141: if (tao->setupcalled) {
142: PetscCall(VecDestroy(&blmP->unprojected_gradient));
143: PetscCall(VecDestroy(&blmP->Xold));
144: PetscCall(VecDestroy(&blmP->Gold));
145: }
146: PetscCall(MatDestroy(&blmP->M));
147: if (blmP->H0) PetscCall(PetscObjectDereference((PetscObject)blmP->H0));
148: PetscCall(PetscFree(tao->data));
149: PetscFunctionReturn(PETSC_SUCCESS);
150: }
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: PetscFunctionBegin;
158: PetscOptionsHeadBegin(PetscOptionsObject, "Limited-memory variable-metric method for bound constrained optimization");
159: PetscCall(PetscOptionsBool("-tao_blmvm_recycle", "enable recycling of the BFGS matrix between subsequent TaoSolve() calls", "", blmP->recycle, &blmP->recycle, NULL));
160: PetscOptionsHeadEnd();
161: PetscCall(MatSetOptionsPrefix(blmP->M, ((PetscObject)tao)->prefix));
162: PetscCall(MatAppendOptionsPrefix(blmP->M, "tao_blmvm_"));
163: PetscCall(MatSetFromOptions(blmP->M));
164: PetscCall(MatIsSPDKnown(blmP->M, &is_set, &is_spd));
165: PetscCheck(is_set && is_spd, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix must be symmetric positive-definite");
166: PetscFunctionReturn(PETSC_SUCCESS);
167: }
169: static PetscErrorCode TaoView_BLMVM(Tao tao, PetscViewer viewer)
170: {
171: TAO_BLMVM *lmP = (TAO_BLMVM *)tao->data;
172: PetscBool isascii;
174: PetscFunctionBegin;
175: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
176: if (isascii) {
177: PetscCall(PetscViewerASCIIPrintf(viewer, "Gradient steps: %" PetscInt_FMT "\n", lmP->grad));
178: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO));
179: PetscCall(MatView(lmP->M, viewer));
180: PetscCall(PetscViewerPopFormat(viewer));
181: }
182: PetscFunctionReturn(PETSC_SUCCESS);
183: }
185: static PetscErrorCode TaoComputeDual_BLMVM(Tao tao, Vec DXL, Vec DXU)
186: {
187: TAO_BLMVM *blm = (TAO_BLMVM *)tao->data;
189: PetscFunctionBegin;
193: PetscCheck(tao->gradient && blm->unprojected_gradient, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Dual variables don't exist yet or no longer exist.");
195: PetscCall(VecCopy(tao->gradient, DXL));
196: PetscCall(VecAXPY(DXL, -1.0, blm->unprojected_gradient));
197: PetscCall(VecSet(DXU, 0.0));
198: PetscCall(VecPointwiseMax(DXL, DXL, DXU));
200: PetscCall(VecCopy(blm->unprojected_gradient, DXU));
201: PetscCall(VecAXPY(DXU, -1.0, tao->gradient));
202: PetscCall(VecAXPY(DXU, 1.0, DXL));
203: PetscFunctionReturn(PETSC_SUCCESS);
204: }
206: /*MC
207: TAOBLMVM - Bounded limited memory variable metric is a quasi-Newton method
208: for nonlinear minimization with bound constraints. It is an extension
209: of `TAOLMVM`
211: Options Database Key:
212: . -tao_lmm_recycle - enable recycling of LMVM information between subsequent `TaoSolve()` calls
214: Level: beginner
216: .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`, `TaoLMVMGetH0()`, `TaoLMVMGetH0KSP()`
217: M*/
218: PETSC_EXTERN PetscErrorCode TaoCreate_BLMVM(Tao tao)
219: {
220: TAO_BLMVM *blmP;
221: const char *morethuente_type = TAOLINESEARCHMT;
223: PetscFunctionBegin;
224: tao->ops->setup = TaoSetup_BLMVM;
225: tao->ops->solve = TaoSolve_BLMVM;
226: tao->ops->view = TaoView_BLMVM;
227: tao->ops->setfromoptions = TaoSetFromOptions_BLMVM;
228: tao->ops->destroy = TaoDestroy_BLMVM;
229: tao->ops->computedual = TaoComputeDual_BLMVM;
231: PetscCall(PetscNew(&blmP));
232: blmP->H0 = NULL;
233: blmP->recycle = PETSC_FALSE;
234: tao->data = (void *)blmP;
236: /* Override default settings (unless already changed) */
237: PetscCall(TaoParametersInitialize(tao));
238: PetscObjectParameterSetDefault(tao, max_it, 2000);
239: PetscObjectParameterSetDefault(tao, max_funcs, 4000);
241: PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
242: PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
243: PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type));
244: PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch, tao));
246: PetscCall(KSPInitializePackage());
247: PetscCall(MatCreate(((PetscObject)tao)->comm, &blmP->M));
248: PetscCall(MatSetType(blmP->M, MATLMVMBFGS));
249: PetscCall(PetscObjectIncrementTabLevel((PetscObject)blmP->M, (PetscObject)tao, 1));
250: PetscFunctionReturn(PETSC_SUCCESS);
251: }
253: /*@
254: TaoLMVMRecycle - Enable/disable recycling of the QN history between subsequent `TaoSolve()` calls.
256: Input Parameters:
257: + tao - the `Tao` solver context
258: - flg - Boolean flag for recycling (`PETSC_TRUE` or `PETSC_FALSE`)
260: Level: intermediate
262: .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`
263: @*/
264: PetscErrorCode TaoLMVMRecycle(Tao tao, PetscBool flg)
265: {
266: TAO_LMVM *lmP;
267: TAO_BLMVM *blmP;
268: PetscBool is_lmvm, is_blmvm;
270: PetscFunctionBegin;
271: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm));
272: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm));
273: if (is_lmvm) {
274: lmP = (TAO_LMVM *)tao->data;
275: lmP->recycle = flg;
276: } else if (is_blmvm) {
277: blmP = (TAO_BLMVM *)tao->data;
278: blmP->recycle = flg;
279: }
280: PetscFunctionReturn(PETSC_SUCCESS);
281: }
283: /*@
284: TaoLMVMSetH0 - Set the initial Hessian for the QN approximation
286: Input Parameters:
287: + tao - the `Tao` solver context
288: - H0 - `Mat` object for the initial Hessian
290: Level: advanced
292: .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`, `TaoLMVMGetH0()`, `TaoLMVMGetH0KSP()`
293: @*/
294: PetscErrorCode TaoLMVMSetH0(Tao tao, Mat H0)
295: {
296: TAO_LMVM *lmP;
297: TAO_BLMVM *blmP;
298: PetscBool is_lmvm, is_blmvm;
300: PetscFunctionBegin;
301: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm));
302: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm));
303: if (is_lmvm) {
304: lmP = (TAO_LMVM *)tao->data;
305: PetscCall(PetscObjectReference((PetscObject)H0));
306: lmP->H0 = H0;
307: } else if (is_blmvm) {
308: blmP = (TAO_BLMVM *)tao->data;
309: PetscCall(PetscObjectReference((PetscObject)H0));
310: blmP->H0 = H0;
311: }
312: PetscFunctionReturn(PETSC_SUCCESS);
313: }
315: /*@
316: TaoLMVMGetH0 - Get the matrix object for the QN initial Hessian
318: Input Parameter:
319: . tao - the `Tao` solver context
321: Output Parameter:
322: . H0 - `Mat` object for the initial Hessian
324: Level: advanced
326: .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`, `TaoLMVMSetH0()`, `TaoLMVMGetH0KSP()`
327: @*/
328: PetscErrorCode TaoLMVMGetH0(Tao tao, Mat *H0)
329: {
330: TAO_LMVM *lmP;
331: TAO_BLMVM *blmP;
332: PetscBool is_lmvm, is_blmvm;
333: Mat M;
335: PetscFunctionBegin;
336: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm));
337: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm));
338: if (is_lmvm) {
339: lmP = (TAO_LMVM *)tao->data;
340: M = lmP->M;
341: } else if (is_blmvm) {
342: blmP = (TAO_BLMVM *)tao->data;
343: M = blmP->M;
344: } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONG, "This routine applies to TAO_LMVM and TAO_BLMVM.");
345: PetscCall(MatLMVMGetJ0(M, H0));
346: PetscFunctionReturn(PETSC_SUCCESS);
347: }
349: /*@
350: TaoLMVMGetH0KSP - Get the iterative solver for applying the inverse of the QN initial Hessian
352: Input Parameter:
353: . tao - the `Tao` solver context
355: Output Parameter:
356: . ksp - `KSP` solver context for the initial Hessian
358: Level: advanced
360: .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`, `TaoLMVMGetH0()`
361: @*/
362: PetscErrorCode TaoLMVMGetH0KSP(Tao tao, KSP *ksp)
363: {
364: TAO_LMVM *lmP;
365: TAO_BLMVM *blmP;
366: PetscBool is_lmvm, is_blmvm;
367: Mat M;
369: PetscFunctionBegin;
370: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm));
371: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm));
372: if (is_lmvm) {
373: lmP = (TAO_LMVM *)tao->data;
374: M = lmP->M;
375: } else if (is_blmvm) {
376: blmP = (TAO_BLMVM *)tao->data;
377: M = blmP->M;
378: } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONG, "This routine applies to TAO_LMVM and TAO_BLMVM.");
379: PetscCall(MatLMVMGetJ0KSP(M, ksp));
380: PetscFunctionReturn(PETSC_SUCCESS);
381: }