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: PetscFunctionBegin;
14: /* Project initial point onto bounds */
15: PetscCall(TaoComputeVariableBounds(tao));
16: PetscCall(VecMedian(tao->XL, tao->solution, tao->XU, tao->solution));
17: PetscCall(TaoLineSearchSetVariableBounds(tao->linesearch, tao->XL, tao->XU));
19: /* Check convergence criteria */
20: PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, blmP->unprojected_gradient));
21: PetscCall(VecBoundGradientProjection(blmP->unprojected_gradient, tao->solution, tao->XL, tao->XU, tao->gradient));
23: PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm));
24: PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Inf or NaN");
26: tao->reason = TAO_CONTINUE_ITERATING;
27: PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
28: PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, stepsize));
29: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
30: if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(PETSC_SUCCESS);
32: /* Set counter for gradient/reset steps */
33: if (!blmP->recycle) {
34: blmP->grad = 0;
35: blmP->reset = 0;
36: PetscCall(MatLMVMReset(blmP->M, PETSC_FALSE));
37: }
39: /* Have not converged; continue with Newton method */
40: while (tao->reason == TAO_CONTINUE_ITERATING) {
41: /* Call general purpose update function */
42: if (tao->ops->update) {
43: PetscUseTypeMethod(tao, update, tao->niter, tao->user_update);
44: PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
45: }
46: /* Compute direction */
47: gnorm2 = gnorm * gnorm;
48: if (gnorm2 == 0.0) gnorm2 = PETSC_MACHINE_EPSILON;
49: if (f == 0.0) {
50: delta = 2.0 / gnorm2;
51: } else {
52: delta = 2.0 * PetscAbsScalar(f) / gnorm2;
53: }
54: PetscCall(MatLMVMSymBroydenSetDelta(blmP->M, delta));
55: PetscCall(MatLMVMUpdate(blmP->M, tao->solution, tao->gradient));
56: PetscCall(MatSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection));
57: PetscCall(VecBoundGradientProjection(tao->stepdirection, tao->solution, tao->XL, tao->XU, tao->gradient));
59: /* Check for success (descent direction) */
60: PetscCall(VecDot(blmP->unprojected_gradient, tao->gradient, &gdx));
61: if (gdx <= 0) {
62: /* Step is not descent or solve was not successful
63: Use steepest descent direction (scaled) */
64: ++blmP->grad;
66: PetscCall(MatLMVMReset(blmP->M, PETSC_FALSE));
67: PetscCall(MatLMVMUpdate(blmP->M, tao->solution, blmP->unprojected_gradient));
68: PetscCall(MatSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection));
69: }
70: PetscCall(VecScale(tao->stepdirection, -1.0));
72: /* Perform the linesearch */
73: fold = f;
74: PetscCall(VecCopy(tao->solution, blmP->Xold));
75: PetscCall(VecCopy(blmP->unprojected_gradient, blmP->Gold));
76: PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0));
77: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, blmP->unprojected_gradient, tao->stepdirection, &stepsize, &ls_status));
78: PetscCall(TaoAddLineSearchCounts(tao));
80: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
81: /* Linesearch failed
82: Reset factors and use scaled (projected) gradient step */
83: ++blmP->reset;
85: f = fold;
86: PetscCall(VecCopy(blmP->Xold, tao->solution));
87: PetscCall(VecCopy(blmP->Gold, blmP->unprojected_gradient));
89: PetscCall(MatLMVMReset(blmP->M, PETSC_FALSE));
90: PetscCall(MatLMVMUpdate(blmP->M, tao->solution, blmP->unprojected_gradient));
91: PetscCall(MatSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection));
92: PetscCall(VecScale(tao->stepdirection, -1.0));
94: /* This may be incorrect; linesearch has values for stepmax and stepmin
95: that should be reset. */
96: PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0));
97: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, blmP->unprojected_gradient, tao->stepdirection, &stepsize, &ls_status));
98: PetscCall(TaoAddLineSearchCounts(tao));
100: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
101: tao->reason = TAO_DIVERGED_LS_FAILURE;
102: break;
103: }
104: }
106: /* Check for converged */
107: PetscCall(VecBoundGradientProjection(blmP->unprojected_gradient, tao->solution, tao->XL, tao->XU, tao->gradient));
108: PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm));
109: PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Not-a-Number");
110: tao->niter++;
111: PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
112: PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, stepsize));
113: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
114: }
115: PetscFunctionReturn(PETSC_SUCCESS);
116: }
118: static PetscErrorCode TaoSetup_BLMVM(Tao tao)
119: {
120: TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data;
122: PetscFunctionBegin;
123: /* Existence of tao->solution checked in TaoSetup() */
124: PetscCall(VecDuplicate(tao->solution, &blmP->Xold));
125: PetscCall(VecDuplicate(tao->solution, &blmP->Gold));
126: PetscCall(VecDuplicate(tao->solution, &blmP->unprojected_gradient));
127: if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
128: if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient));
129: /* Allocate matrix for the limited memory approximation */
130: PetscCall(MatLMVMAllocate(blmP->M, tao->solution, blmP->unprojected_gradient));
132: /* If the user has set a matrix to solve as the initial H0, set the options prefix here, and set up the KSP */
133: if (blmP->H0) PetscCall(MatLMVMSetJ0(blmP->M, blmP->H0));
134: PetscFunctionReturn(PETSC_SUCCESS);
135: }
137: /* ---------------------------------------------------------- */
138: static PetscErrorCode TaoDestroy_BLMVM(Tao tao)
139: {
140: TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data;
142: PetscFunctionBegin;
143: if (tao->setupcalled) {
144: PetscCall(VecDestroy(&blmP->unprojected_gradient));
145: PetscCall(VecDestroy(&blmP->Xold));
146: PetscCall(VecDestroy(&blmP->Gold));
147: }
148: PetscCall(MatDestroy(&blmP->M));
149: if (blmP->H0) PetscCall(PetscObjectDereference((PetscObject)blmP->H0));
150: PetscCall(PetscFree(tao->data));
151: PetscFunctionReturn(PETSC_SUCCESS);
152: }
154: /*------------------------------------------------------------*/
155: static PetscErrorCode TaoSetFromOptions_BLMVM(Tao tao, PetscOptionItems *PetscOptionsObject)
156: {
157: TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data;
158: PetscBool is_spd, is_set;
160: PetscFunctionBegin;
161: PetscOptionsHeadBegin(PetscOptionsObject, "Limited-memory variable-metric method for bound constrained optimization");
162: PetscCall(PetscOptionsBool("-tao_blmvm_recycle", "enable recycling of the BFGS matrix between subsequent TaoSolve() calls", "", blmP->recycle, &blmP->recycle, NULL));
163: PetscOptionsHeadEnd();
164: PetscCall(MatSetOptionsPrefix(blmP->M, ((PetscObject)tao)->prefix));
165: PetscCall(MatAppendOptionsPrefix(blmP->M, "tao_blmvm_"));
166: PetscCall(MatSetFromOptions(blmP->M));
167: PetscCall(MatIsSPDKnown(blmP->M, &is_set, &is_spd));
168: PetscCheck(is_set && is_spd, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix must be symmetric positive-definite");
169: PetscFunctionReturn(PETSC_SUCCESS);
170: }
172: /*------------------------------------------------------------*/
173: static PetscErrorCode TaoView_BLMVM(Tao tao, PetscViewer viewer)
174: {
175: TAO_BLMVM *lmP = (TAO_BLMVM *)tao->data;
176: PetscBool isascii;
178: PetscFunctionBegin;
179: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
180: if (isascii) {
181: PetscCall(PetscViewerASCIIPrintf(viewer, "Gradient steps: %" PetscInt_FMT "\n", lmP->grad));
182: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO));
183: PetscCall(MatView(lmP->M, viewer));
184: PetscCall(PetscViewerPopFormat(viewer));
185: }
186: PetscFunctionReturn(PETSC_SUCCESS);
187: }
189: static PetscErrorCode TaoComputeDual_BLMVM(Tao tao, Vec DXL, Vec DXU)
190: {
191: TAO_BLMVM *blm = (TAO_BLMVM *)tao->data;
193: PetscFunctionBegin;
197: PetscCheck(tao->gradient && blm->unprojected_gradient, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Dual variables don't exist yet or no longer exist.");
199: PetscCall(VecCopy(tao->gradient, DXL));
200: PetscCall(VecAXPY(DXL, -1.0, blm->unprojected_gradient));
201: PetscCall(VecSet(DXU, 0.0));
202: PetscCall(VecPointwiseMax(DXL, DXL, DXU));
204: PetscCall(VecCopy(blm->unprojected_gradient, DXU));
205: PetscCall(VecAXPY(DXU, -1.0, tao->gradient));
206: PetscCall(VecAXPY(DXU, 1.0, DXL));
207: PetscFunctionReturn(PETSC_SUCCESS);
208: }
210: /* ---------------------------------------------------------- */
211: /*MC
212: TAOBLMVM - Bounded limited memory variable metric is a quasi-Newton method
213: for nonlinear minimization with bound constraints. It is an extension
214: of `TAOLMVM`
216: Options Database Key:
217: . -tao_lmm_recycle - enable recycling of LMVM information between subsequent `TaoSolve()` calls
219: Level: beginner
221: .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`, `TaoLMVMGetH0()`, `TaoLMVMGetH0KSP()`
222: M*/
223: PETSC_EXTERN PetscErrorCode TaoCreate_BLMVM(Tao tao)
224: {
225: TAO_BLMVM *blmP;
226: const char *morethuente_type = TAOLINESEARCHMT;
228: PetscFunctionBegin;
229: tao->ops->setup = TaoSetup_BLMVM;
230: tao->ops->solve = TaoSolve_BLMVM;
231: tao->ops->view = TaoView_BLMVM;
232: tao->ops->setfromoptions = TaoSetFromOptions_BLMVM;
233: tao->ops->destroy = TaoDestroy_BLMVM;
234: tao->ops->computedual = TaoComputeDual_BLMVM;
236: PetscCall(PetscNew(&blmP));
237: blmP->H0 = NULL;
238: blmP->recycle = PETSC_FALSE;
239: tao->data = (void *)blmP;
241: /* Override default settings (unless already changed) */
242: if (!tao->max_it_changed) tao->max_it = 2000;
243: if (!tao->max_funcs_changed) tao->max_funcs = 4000;
245: PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
246: PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
247: PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type));
248: PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch, tao));
250: PetscCall(KSPInitializePackage());
251: PetscCall(MatCreate(((PetscObject)tao)->comm, &blmP->M));
252: PetscCall(MatSetType(blmP->M, MATLMVMBFGS));
253: PetscCall(PetscObjectIncrementTabLevel((PetscObject)blmP->M, (PetscObject)tao, 1));
254: PetscFunctionReturn(PETSC_SUCCESS);
255: }
257: /*@
258: TaoLMVMRecycle - Enable/disable recycling of the QN history between subsequent `TaoSolve()` calls.
260: Input Parameters:
261: + tao - the `Tao` solver context
262: - flg - Boolean flag for recycling (`PETSC_TRUE` or `PETSC_FALSE`)
264: Level: intermediate
266: .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`
267: @*/
268: PetscErrorCode TaoLMVMRecycle(Tao tao, PetscBool flg)
269: {
270: TAO_LMVM *lmP;
271: TAO_BLMVM *blmP;
272: PetscBool is_lmvm, is_blmvm;
274: PetscFunctionBegin;
275: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm));
276: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm));
277: if (is_lmvm) {
278: lmP = (TAO_LMVM *)tao->data;
279: lmP->recycle = flg;
280: } else if (is_blmvm) {
281: blmP = (TAO_BLMVM *)tao->data;
282: blmP->recycle = flg;
283: }
284: PetscFunctionReturn(PETSC_SUCCESS);
285: }
287: /*@
288: TaoLMVMSetH0 - Set the initial Hessian for the QN approximation
290: Input Parameters:
291: + tao - the `Tao` solver context
292: - H0 - `Mat` object for the initial Hessian
294: Level: advanced
296: .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`, `TaoLMVMGetH0()`, `TaoLMVMGetH0KSP()`
297: @*/
298: PetscErrorCode TaoLMVMSetH0(Tao tao, Mat H0)
299: {
300: TAO_LMVM *lmP;
301: TAO_BLMVM *blmP;
302: PetscBool is_lmvm, is_blmvm;
304: PetscFunctionBegin;
305: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm));
306: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm));
307: if (is_lmvm) {
308: lmP = (TAO_LMVM *)tao->data;
309: PetscCall(PetscObjectReference((PetscObject)H0));
310: lmP->H0 = H0;
311: } else if (is_blmvm) {
312: blmP = (TAO_BLMVM *)tao->data;
313: PetscCall(PetscObjectReference((PetscObject)H0));
314: blmP->H0 = H0;
315: }
316: PetscFunctionReturn(PETSC_SUCCESS);
317: }
319: /*@
320: TaoLMVMGetH0 - Get the matrix object for the QN initial Hessian
322: Input Parameter:
323: . tao - the `Tao` solver context
325: Output Parameter:
326: . H0 - `Mat` object for the initial Hessian
328: Level: advanced
330: .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`, `TaoLMVMSetH0()`, `TaoLMVMGetH0KSP()`
331: @*/
332: PetscErrorCode TaoLMVMGetH0(Tao tao, Mat *H0)
333: {
334: TAO_LMVM *lmP;
335: TAO_BLMVM *blmP;
336: PetscBool is_lmvm, is_blmvm;
337: Mat M;
339: PetscFunctionBegin;
340: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm));
341: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm));
342: if (is_lmvm) {
343: lmP = (TAO_LMVM *)tao->data;
344: M = lmP->M;
345: } else if (is_blmvm) {
346: blmP = (TAO_BLMVM *)tao->data;
347: M = blmP->M;
348: } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONG, "This routine applies to TAO_LMVM and TAO_BLMVM.");
349: PetscCall(MatLMVMGetJ0(M, H0));
350: PetscFunctionReturn(PETSC_SUCCESS);
351: }
353: /*@
354: TaoLMVMGetH0KSP - Get the iterative solver for applying the inverse of the QN initial Hessian
356: Input Parameter:
357: . tao - the `Tao` solver context
359: Output Parameter:
360: . ksp - `KSP` solver context for the initial Hessian
362: Level: advanced
364: .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`, `TaoLMVMGetH0()`
365: @*/
366: PetscErrorCode TaoLMVMGetH0KSP(Tao tao, KSP *ksp)
367: {
368: TAO_LMVM *lmP;
369: TAO_BLMVM *blmP;
370: PetscBool is_lmvm, is_blmvm;
371: Mat M;
373: PetscFunctionBegin;
374: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm));
375: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm));
376: if (is_lmvm) {
377: lmP = (TAO_LMVM *)tao->data;
378: M = lmP->M;
379: } else if (is_blmvm) {
380: blmP = (TAO_BLMVM *)tao->data;
381: M = blmP->M;
382: } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONG, "This routine applies to TAO_LMVM and TAO_BLMVM.");
383: PetscCall(MatLMVMGetJ0KSP(M, ksp));
384: PetscFunctionReturn(PETSC_SUCCESS);
385: }