1 #include <petsctaolinesearch.h> /*I "petsctaolinesearch.h" I*/
2 #include <../src/tao/unconstrained/impls/lmvm/lmvm.h>
3 #include <../src/tao/bound/impls/blmvm/blmvm.h>
4
TaoSolve_BLMVM(Tao tao)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;
11
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));
17
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));
21
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");
24
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);
30
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 }
37
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));
57
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;
64
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));
70
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));
78
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;
83
84 f = fold;
85 PetscCall(VecCopy(blmP->Xold, tao->solution));
86 PetscCall(VecCopy(blmP->Gold, blmP->unprojected_gradient));
87
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));
92
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));
98
99 if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
100 tao->reason = TAO_DIVERGED_LS_FAILURE;
101 break;
102 }
103 }
104
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 }
116
TaoSetup_BLMVM(Tao tao)117 static PetscErrorCode TaoSetup_BLMVM(Tao tao)
118 {
119 TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data;
120
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));
130
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 }
135
TaoDestroy_BLMVM(Tao tao)136 static PetscErrorCode TaoDestroy_BLMVM(Tao tao)
137 {
138 TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data;
139
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 }
151
TaoSetFromOptions_BLMVM(Tao tao,PetscOptionItems PetscOptionsObject)152 static PetscErrorCode TaoSetFromOptions_BLMVM(Tao tao, PetscOptionItems PetscOptionsObject)
153 {
154 TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data;
155 PetscBool is_spd, is_set;
156
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 }
168
TaoView_BLMVM(Tao tao,PetscViewer viewer)169 static PetscErrorCode TaoView_BLMVM(Tao tao, PetscViewer viewer)
170 {
171 TAO_BLMVM *lmP = (TAO_BLMVM *)tao->data;
172 PetscBool isascii;
173
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 }
184
TaoComputeDual_BLMVM(Tao tao,Vec DXL,Vec DXU)185 static PetscErrorCode TaoComputeDual_BLMVM(Tao tao, Vec DXL, Vec DXU)
186 {
187 TAO_BLMVM *blm = (TAO_BLMVM *)tao->data;
188
189 PetscFunctionBegin;
190 PetscValidHeaderSpecific(tao, TAO_CLASSID, 1);
191 PetscValidHeaderSpecific(DXL, VEC_CLASSID, 2);
192 PetscValidHeaderSpecific(DXU, VEC_CLASSID, 3);
193 PetscCheck(tao->gradient && blm->unprojected_gradient, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Dual variables don't exist yet or no longer exist.");
194
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));
199
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 }
205
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`
210
211 Options Database Key:
212 . -tao_lmm_recycle - enable recycling of LMVM information between subsequent `TaoSolve()` calls
213
214 Level: beginner
215
216 .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`, `TaoLMVMGetH0()`, `TaoLMVMGetH0KSP()`
217 M*/
TaoCreate_BLMVM(Tao tao)218 PETSC_EXTERN PetscErrorCode TaoCreate_BLMVM(Tao tao)
219 {
220 TAO_BLMVM *blmP;
221 const char *morethuente_type = TAOLINESEARCHMT;
222
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;
230
231 PetscCall(PetscNew(&blmP));
232 blmP->H0 = NULL;
233 blmP->recycle = PETSC_FALSE;
234 tao->data = (void *)blmP;
235
236 /* Override default settings (unless already changed) */
237 PetscCall(TaoParametersInitialize(tao));
238 PetscObjectParameterSetDefault(tao, max_it, 2000);
239 PetscObjectParameterSetDefault(tao, max_funcs, 4000);
240
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));
245
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 }
252
253 /*@
254 TaoLMVMRecycle - Enable/disable recycling of the QN history between subsequent `TaoSolve()` calls.
255
256 Input Parameters:
257 + tao - the `Tao` solver context
258 - flg - Boolean flag for recycling (`PETSC_TRUE` or `PETSC_FALSE`)
259
260 Level: intermediate
261
262 .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`
263 @*/
TaoLMVMRecycle(Tao tao,PetscBool flg)264 PetscErrorCode TaoLMVMRecycle(Tao tao, PetscBool flg)
265 {
266 TAO_LMVM *lmP;
267 TAO_BLMVM *blmP;
268 PetscBool is_lmvm, is_blmvm;
269
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 }
282
283 /*@
284 TaoLMVMSetH0 - Set the initial Hessian for the QN approximation
285
286 Input Parameters:
287 + tao - the `Tao` solver context
288 - H0 - `Mat` object for the initial Hessian
289
290 Level: advanced
291
292 .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`, `TaoLMVMGetH0()`, `TaoLMVMGetH0KSP()`
293 @*/
TaoLMVMSetH0(Tao tao,Mat H0)294 PetscErrorCode TaoLMVMSetH0(Tao tao, Mat H0)
295 {
296 TAO_LMVM *lmP;
297 TAO_BLMVM *blmP;
298 PetscBool is_lmvm, is_blmvm;
299
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 }
314
315 /*@
316 TaoLMVMGetH0 - Get the matrix object for the QN initial Hessian
317
318 Input Parameter:
319 . tao - the `Tao` solver context
320
321 Output Parameter:
322 . H0 - `Mat` object for the initial Hessian
323
324 Level: advanced
325
326 .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`, `TaoLMVMSetH0()`, `TaoLMVMGetH0KSP()`
327 @*/
TaoLMVMGetH0(Tao tao,Mat * H0)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;
334
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 }
348
349 /*@
350 TaoLMVMGetH0KSP - Get the iterative solver for applying the inverse of the QN initial Hessian
351
352 Input Parameter:
353 . tao - the `Tao` solver context
354
355 Output Parameter:
356 . ksp - `KSP` solver context for the initial Hessian
357
358 Level: advanced
359
360 .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`, `TaoLMVMGetH0()`
361 @*/
TaoLMVMGetH0KSP(Tao tao,KSP * ksp)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;
368
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 }
382