xref: /petsc/src/tao/unconstrained/impls/ntl/ntl.c (revision d2522c19e8fa9bca20aaca277941d9a63e71db6a) 
1 #include <../src/tao/unconstrained/impls/ntl/ntlimpl.h>
2 
3 #include <petscksp.h>
4 
5 #define NTL_INIT_CONSTANT      0
6 #define NTL_INIT_DIRECTION     1
7 #define NTL_INIT_INTERPOLATION 2
8 #define NTL_INIT_TYPES         3
9 
10 #define NTL_UPDATE_REDUCTION     0
11 #define NTL_UPDATE_INTERPOLATION 1
12 #define NTL_UPDATE_TYPES         2
13 
14 static const char *NTL_INIT[64] = {"constant", "direction", "interpolation"};
15 
16 static const char *NTL_UPDATE[64] = {"reduction", "interpolation"};
17 
18 /* Implements Newton's Method with a trust-region, line-search approach for
19    solving unconstrained minimization problems.  A More'-Thuente line search
20    is used to guarantee that the bfgs preconditioner remains positive
21    definite. */
22 
23 #define NTL_NEWTON          0
24 #define NTL_BFGS            1
25 #define NTL_SCALED_GRADIENT 2
26 #define NTL_GRADIENT        3
27 
28 static PetscErrorCode TaoSolve_NTL(Tao tao) {
29   TAO_NTL                     *tl = (TAO_NTL *)tao->data;
30   KSPType                      ksp_type;
31   PetscBool                    is_nash, is_stcg, is_gltr, is_bfgs, is_jacobi, is_symmetric, sym_set;
32   KSPConvergedReason           ksp_reason;
33   PC                           pc;
34   TaoLineSearchConvergedReason ls_reason;
35 
36   PetscReal fmin, ftrial, prered, actred, kappa, sigma;
37   PetscReal tau, tau_1, tau_2, tau_max, tau_min, max_radius;
38   PetscReal f, fold, gdx, gnorm;
39   PetscReal step = 1.0;
40 
41   PetscReal norm_d = 0.0;
42   PetscInt  stepType;
43   PetscInt  its;
44 
45   PetscInt bfgsUpdates = 0;
46   PetscInt needH;
47 
48   PetscInt i_max = 5;
49   PetscInt j_max = 1;
50   PetscInt i, j, n, N;
51 
52   PetscInt tr_reject;
53 
54   PetscFunctionBegin;
55   if (tao->XL || tao->XU || tao->ops->computebounds) { PetscCall(PetscInfo(tao, "WARNING: Variable bounds have been set but will be ignored by ntl algorithm\n")); }
56 
57   PetscCall(KSPGetType(tao->ksp, &ksp_type));
58   PetscCall(PetscStrcmp(ksp_type, KSPNASH, &is_nash));
59   PetscCall(PetscStrcmp(ksp_type, KSPSTCG, &is_stcg));
60   PetscCall(PetscStrcmp(ksp_type, KSPGLTR, &is_gltr));
61   PetscCheck(is_nash || is_stcg || is_gltr, PetscObjectComm((PetscObject)tao), PETSC_ERR_SUP, "TAO_NTR requires nash, stcg, or gltr for the KSP");
62 
63   /* Initialize the radius and modify if it is too large or small */
64   tao->trust = tao->trust0;
65   tao->trust = PetscMax(tao->trust, tl->min_radius);
66   tao->trust = PetscMin(tao->trust, tl->max_radius);
67 
68   /* Allocate the vectors needed for the BFGS approximation */
69   PetscCall(KSPGetPC(tao->ksp, &pc));
70   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &is_bfgs));
71   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &is_jacobi));
72   if (is_bfgs) {
73     tl->bfgs_pre = pc;
74     PetscCall(PCLMVMGetMatLMVM(tl->bfgs_pre, &tl->M));
75     PetscCall(VecGetLocalSize(tao->solution, &n));
76     PetscCall(VecGetSize(tao->solution, &N));
77     PetscCall(MatSetSizes(tl->M, n, n, N, N));
78     PetscCall(MatLMVMAllocate(tl->M, tao->solution, tao->gradient));
79     PetscCall(MatIsSymmetricKnown(tl->M, &sym_set, &is_symmetric));
80     PetscCheck(sym_set && is_symmetric, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix in the LMVM preconditioner must be symmetric.");
81   } else if (is_jacobi) PetscCall(PCJacobiSetUseAbs(pc, PETSC_TRUE));
82 
83   /* Check convergence criteria */
84   PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
85   PetscCall(VecNorm(tao->gradient, NORM_2, &gnorm));
86   PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Inf or NaN");
87   needH = 1;
88 
89   tao->reason = TAO_CONTINUE_ITERATING;
90   PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
91   PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step));
92   PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
93   if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0);
94 
95   /* Initialize trust-region radius */
96   switch (tl->init_type) {
97   case NTL_INIT_CONSTANT:
98     /* Use the initial radius specified */
99     break;
100 
101   case NTL_INIT_INTERPOLATION:
102     /* Use the initial radius specified */
103     max_radius = 0.0;
104 
105     for (j = 0; j < j_max; ++j) {
106       fmin  = f;
107       sigma = 0.0;
108 
109       if (needH) {
110         PetscCall(TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre));
111         needH = 0;
112       }
113 
114       for (i = 0; i < i_max; ++i) {
115         PetscCall(VecCopy(tao->solution, tl->W));
116         PetscCall(VecAXPY(tl->W, -tao->trust / gnorm, tao->gradient));
117 
118         PetscCall(TaoComputeObjective(tao, tl->W, &ftrial));
119         if (PetscIsInfOrNanReal(ftrial)) {
120           tau = tl->gamma1_i;
121         } else {
122           if (ftrial < fmin) {
123             fmin  = ftrial;
124             sigma = -tao->trust / gnorm;
125           }
126 
127           PetscCall(MatMult(tao->hessian, tao->gradient, tao->stepdirection));
128           PetscCall(VecDot(tao->gradient, tao->stepdirection, &prered));
129 
130           prered = tao->trust * (gnorm - 0.5 * tao->trust * prered / (gnorm * gnorm));
131           actred = f - ftrial;
132           if ((PetscAbsScalar(actred) <= tl->epsilon) && (PetscAbsScalar(prered) <= tl->epsilon)) {
133             kappa = 1.0;
134           } else {
135             kappa = actred / prered;
136           }
137 
138           tau_1   = tl->theta_i * gnorm * tao->trust / (tl->theta_i * gnorm * tao->trust + (1.0 - tl->theta_i) * prered - actred);
139           tau_2   = tl->theta_i * gnorm * tao->trust / (tl->theta_i * gnorm * tao->trust - (1.0 + tl->theta_i) * prered + actred);
140           tau_min = PetscMin(tau_1, tau_2);
141           tau_max = PetscMax(tau_1, tau_2);
142 
143           if (PetscAbsScalar(kappa - (PetscReal)1.0) <= tl->mu1_i) {
144             /* Great agreement */
145             max_radius = PetscMax(max_radius, tao->trust);
146 
147             if (tau_max < 1.0) {
148               tau = tl->gamma3_i;
149             } else if (tau_max > tl->gamma4_i) {
150               tau = tl->gamma4_i;
151             } else if (tau_1 >= 1.0 && tau_1 <= tl->gamma4_i && tau_2 < 1.0) {
152               tau = tau_1;
153             } else if (tau_2 >= 1.0 && tau_2 <= tl->gamma4_i && tau_1 < 1.0) {
154               tau = tau_2;
155             } else {
156               tau = tau_max;
157             }
158           } else if (PetscAbsScalar(kappa - (PetscReal)1.0) <= tl->mu2_i) {
159             /* Good agreement */
160             max_radius = PetscMax(max_radius, tao->trust);
161 
162             if (tau_max < tl->gamma2_i) {
163               tau = tl->gamma2_i;
164             } else if (tau_max > tl->gamma3_i) {
165               tau = tl->gamma3_i;
166             } else {
167               tau = tau_max;
168             }
169           } else {
170             /* Not good agreement */
171             if (tau_min > 1.0) {
172               tau = tl->gamma2_i;
173             } else if (tau_max < tl->gamma1_i) {
174               tau = tl->gamma1_i;
175             } else if ((tau_min < tl->gamma1_i) && (tau_max >= 1.0)) {
176               tau = tl->gamma1_i;
177             } else if ((tau_1 >= tl->gamma1_i) && (tau_1 < 1.0) && ((tau_2 < tl->gamma1_i) || (tau_2 >= 1.0))) {
178               tau = tau_1;
179             } else if ((tau_2 >= tl->gamma1_i) && (tau_2 < 1.0) && ((tau_1 < tl->gamma1_i) || (tau_2 >= 1.0))) {
180               tau = tau_2;
181             } else {
182               tau = tau_max;
183             }
184           }
185         }
186         tao->trust = tau * tao->trust;
187       }
188 
189       if (fmin < f) {
190         f = fmin;
191         PetscCall(VecAXPY(tao->solution, sigma, tao->gradient));
192         PetscCall(TaoComputeGradient(tao, tao->solution, tao->gradient));
193 
194         PetscCall(VecNorm(tao->gradient, NORM_2, &gnorm));
195         PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Inf or NaN");
196         needH = 1;
197 
198         PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
199         PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step));
200         PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
201         if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0);
202       }
203     }
204     tao->trust = PetscMax(tao->trust, max_radius);
205 
206     /* Modify the radius if it is too large or small */
207     tao->trust = PetscMax(tao->trust, tl->min_radius);
208     tao->trust = PetscMin(tao->trust, tl->max_radius);
209     break;
210 
211   default:
212     /* Norm of the first direction will initialize radius */
213     tao->trust = 0.0;
214     break;
215   }
216 
217   /* Set counter for gradient/reset steps */
218   tl->ntrust = 0;
219   tl->newt   = 0;
220   tl->bfgs   = 0;
221   tl->grad   = 0;
222 
223   /* Have not converged; continue with Newton method */
224   while (tao->reason == TAO_CONTINUE_ITERATING) {
225     /* Call general purpose update function */
226     PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
227     ++tao->niter;
228     tao->ksp_its = 0;
229     /* Compute the Hessian */
230     if (needH) PetscCall(TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre));
231 
232     if (tl->bfgs_pre) {
233       /* Update the limited memory preconditioner */
234       PetscCall(MatLMVMUpdate(tl->M, tao->solution, tao->gradient));
235       ++bfgsUpdates;
236     }
237     PetscCall(KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre));
238     /* Solve the Newton system of equations */
239     PetscCall(KSPCGSetRadius(tao->ksp, tl->max_radius));
240     PetscCall(KSPSolve(tao->ksp, tao->gradient, tao->stepdirection));
241     PetscCall(KSPGetIterationNumber(tao->ksp, &its));
242     tao->ksp_its += its;
243     tao->ksp_tot_its += its;
244     PetscCall(KSPCGGetNormD(tao->ksp, &norm_d));
245 
246     if (0.0 == tao->trust) {
247       /* Radius was uninitialized; use the norm of the direction */
248       if (norm_d > 0.0) {
249         tao->trust = norm_d;
250 
251         /* Modify the radius if it is too large or small */
252         tao->trust = PetscMax(tao->trust, tl->min_radius);
253         tao->trust = PetscMin(tao->trust, tl->max_radius);
254       } else {
255         /* The direction was bad; set radius to default value and re-solve
256            the trust-region subproblem to get a direction */
257         tao->trust = tao->trust0;
258 
259         /* Modify the radius if it is too large or small */
260         tao->trust = PetscMax(tao->trust, tl->min_radius);
261         tao->trust = PetscMin(tao->trust, tl->max_radius);
262 
263         PetscCall(KSPCGSetRadius(tao->ksp, tl->max_radius));
264         PetscCall(KSPSolve(tao->ksp, tao->gradient, tao->stepdirection));
265         PetscCall(KSPGetIterationNumber(tao->ksp, &its));
266         tao->ksp_its += its;
267         tao->ksp_tot_its += its;
268         PetscCall(KSPCGGetNormD(tao->ksp, &norm_d));
269 
270         PetscCheck(norm_d != 0.0, PetscObjectComm((PetscObject)tao), PETSC_ERR_PLIB, "Initial direction zero");
271       }
272     }
273 
274     PetscCall(VecScale(tao->stepdirection, -1.0));
275     PetscCall(KSPGetConvergedReason(tao->ksp, &ksp_reason));
276     if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) && (tl->bfgs_pre)) {
277       /* Preconditioner is numerically indefinite; reset the
278          approximate if using BFGS preconditioning. */
279       PetscCall(MatLMVMReset(tl->M, PETSC_FALSE));
280       PetscCall(MatLMVMUpdate(tl->M, tao->solution, tao->gradient));
281       bfgsUpdates = 1;
282     }
283 
284     /* Check trust-region reduction conditions */
285     tr_reject = 0;
286     if (NTL_UPDATE_REDUCTION == tl->update_type) {
287       /* Get predicted reduction */
288       PetscCall(KSPCGGetObjFcn(tao->ksp, &prered));
289       if (prered >= 0.0) {
290         /* The predicted reduction has the wrong sign.  This cannot
291            happen in infinite precision arithmetic.  Step should
292            be rejected! */
293         tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d);
294         tr_reject  = 1;
295       } else {
296         /* Compute trial step and function value */
297         PetscCall(VecCopy(tao->solution, tl->W));
298         PetscCall(VecAXPY(tl->W, 1.0, tao->stepdirection));
299         PetscCall(TaoComputeObjective(tao, tl->W, &ftrial));
300 
301         if (PetscIsInfOrNanReal(ftrial)) {
302           tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d);
303           tr_reject  = 1;
304         } else {
305           /* Compute and actual reduction */
306           actred = f - ftrial;
307           prered = -prered;
308           if ((PetscAbsScalar(actred) <= tl->epsilon) && (PetscAbsScalar(prered) <= tl->epsilon)) {
309             kappa = 1.0;
310           } else {
311             kappa = actred / prered;
312           }
313 
314           /* Accept of reject the step and update radius */
315           if (kappa < tl->eta1) {
316             /* Reject the step */
317             tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d);
318             tr_reject  = 1;
319           } else {
320             /* Accept the step */
321             if (kappa < tl->eta2) {
322               /* Marginal bad step */
323               tao->trust = tl->alpha2 * PetscMin(tao->trust, norm_d);
324             } else if (kappa < tl->eta3) {
325               /* Reasonable step */
326               tao->trust = tl->alpha3 * tao->trust;
327             } else if (kappa < tl->eta4) {
328               /* Good step */
329               tao->trust = PetscMax(tl->alpha4 * norm_d, tao->trust);
330             } else {
331               /* Very good step */
332               tao->trust = PetscMax(tl->alpha5 * norm_d, tao->trust);
333             }
334           }
335         }
336       }
337     } else {
338       /* Get predicted reduction */
339       PetscCall(KSPCGGetObjFcn(tao->ksp, &prered));
340       if (prered >= 0.0) {
341         /* The predicted reduction has the wrong sign.  This cannot
342            happen in infinite precision arithmetic.  Step should
343            be rejected! */
344         tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
345         tr_reject  = 1;
346       } else {
347         PetscCall(VecCopy(tao->solution, tl->W));
348         PetscCall(VecAXPY(tl->W, 1.0, tao->stepdirection));
349         PetscCall(TaoComputeObjective(tao, tl->W, &ftrial));
350         if (PetscIsInfOrNanReal(ftrial)) {
351           tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
352           tr_reject  = 1;
353         } else {
354           PetscCall(VecDot(tao->gradient, tao->stepdirection, &gdx));
355 
356           actred = f - ftrial;
357           prered = -prered;
358           if ((PetscAbsScalar(actred) <= tl->epsilon) && (PetscAbsScalar(prered) <= tl->epsilon)) {
359             kappa = 1.0;
360           } else {
361             kappa = actred / prered;
362           }
363 
364           tau_1   = tl->theta * gdx / (tl->theta * gdx - (1.0 - tl->theta) * prered + actred);
365           tau_2   = tl->theta * gdx / (tl->theta * gdx + (1.0 + tl->theta) * prered - actred);
366           tau_min = PetscMin(tau_1, tau_2);
367           tau_max = PetscMax(tau_1, tau_2);
368 
369           if (kappa >= 1.0 - tl->mu1) {
370             /* Great agreement; accept step and update radius */
371             if (tau_max < 1.0) {
372               tao->trust = PetscMax(tao->trust, tl->gamma3 * norm_d);
373             } else if (tau_max > tl->gamma4) {
374               tao->trust = PetscMax(tao->trust, tl->gamma4 * norm_d);
375             } else {
376               tao->trust = PetscMax(tao->trust, tau_max * norm_d);
377             }
378           } else if (kappa >= 1.0 - tl->mu2) {
379             /* Good agreement */
380 
381             if (tau_max < tl->gamma2) {
382               tao->trust = tl->gamma2 * PetscMin(tao->trust, norm_d);
383             } else if (tau_max > tl->gamma3) {
384               tao->trust = PetscMax(tao->trust, tl->gamma3 * norm_d);
385             } else if (tau_max < 1.0) {
386               tao->trust = tau_max * PetscMin(tao->trust, norm_d);
387             } else {
388               tao->trust = PetscMax(tao->trust, tau_max * norm_d);
389             }
390           } else {
391             /* Not good agreement */
392             if (tau_min > 1.0) {
393               tao->trust = tl->gamma2 * PetscMin(tao->trust, norm_d);
394             } else if (tau_max < tl->gamma1) {
395               tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
396             } else if ((tau_min < tl->gamma1) && (tau_max >= 1.0)) {
397               tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
398             } else if ((tau_1 >= tl->gamma1) && (tau_1 < 1.0) && ((tau_2 < tl->gamma1) || (tau_2 >= 1.0))) {
399               tao->trust = tau_1 * PetscMin(tao->trust, norm_d);
400             } else if ((tau_2 >= tl->gamma1) && (tau_2 < 1.0) && ((tau_1 < tl->gamma1) || (tau_2 >= 1.0))) {
401               tao->trust = tau_2 * PetscMin(tao->trust, norm_d);
402             } else {
403               tao->trust = tau_max * PetscMin(tao->trust, norm_d);
404             }
405             tr_reject = 1;
406           }
407         }
408       }
409     }
410 
411     if (tr_reject) {
412       /* The trust-region constraints rejected the step.  Apply a linesearch.
413          Check for descent direction. */
414       PetscCall(VecDot(tao->stepdirection, tao->gradient, &gdx));
415       if ((gdx >= 0.0) || PetscIsInfOrNanReal(gdx)) {
416         /* Newton step is not descent or direction produced Inf or NaN */
417 
418         if (!tl->bfgs_pre) {
419           /* We don't have the bfgs matrix around and updated
420              Must use gradient direction in this case */
421           PetscCall(VecCopy(tao->gradient, tao->stepdirection));
422           PetscCall(VecScale(tao->stepdirection, -1.0));
423           ++tl->grad;
424           stepType = NTL_GRADIENT;
425         } else {
426           /* Attempt to use the BFGS direction */
427           PetscCall(MatSolve(tl->M, tao->gradient, tao->stepdirection));
428           PetscCall(VecScale(tao->stepdirection, -1.0));
429 
430           /* Check for success (descent direction) */
431           PetscCall(VecDot(tao->stepdirection, tao->gradient, &gdx));
432           if ((gdx >= 0) || PetscIsInfOrNanReal(gdx)) {
433             /* BFGS direction is not descent or direction produced not a number
434                We can assert bfgsUpdates > 1 in this case because
435                the first solve produces the scaled gradient direction,
436                which is guaranteed to be descent */
437 
438             /* Use steepest descent direction (scaled) */
439             PetscCall(MatLMVMReset(tl->M, PETSC_FALSE));
440             PetscCall(MatLMVMUpdate(tl->M, tao->solution, tao->gradient));
441             PetscCall(MatSolve(tl->M, tao->gradient, tao->stepdirection));
442             PetscCall(VecScale(tao->stepdirection, -1.0));
443 
444             bfgsUpdates = 1;
445             ++tl->grad;
446             stepType = NTL_GRADIENT;
447           } else {
448             PetscCall(MatLMVMGetUpdateCount(tl->M, &bfgsUpdates));
449             if (1 == bfgsUpdates) {
450               /* The first BFGS direction is always the scaled gradient */
451               ++tl->grad;
452               stepType = NTL_GRADIENT;
453             } else {
454               ++tl->bfgs;
455               stepType = NTL_BFGS;
456             }
457           }
458         }
459       } else {
460         /* Computed Newton step is descent */
461         ++tl->newt;
462         stepType = NTL_NEWTON;
463       }
464 
465       /* Perform the linesearch */
466       fold = f;
467       PetscCall(VecCopy(tao->solution, tl->Xold));
468       PetscCall(VecCopy(tao->gradient, tl->Gold));
469 
470       step = 1.0;
471       PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_reason));
472       PetscCall(TaoAddLineSearchCounts(tao));
473 
474       while (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER && stepType != NTL_GRADIENT) { /* Linesearch failed */
475         /* Linesearch failed */
476         f = fold;
477         PetscCall(VecCopy(tl->Xold, tao->solution));
478         PetscCall(VecCopy(tl->Gold, tao->gradient));
479 
480         switch (stepType) {
481         case NTL_NEWTON:
482           /* Failed to obtain acceptable iterate with Newton step */
483 
484           if (tl->bfgs_pre) {
485             /* We don't have the bfgs matrix around and being updated
486                Must use gradient direction in this case */
487             PetscCall(VecCopy(tao->gradient, tao->stepdirection));
488             ++tl->grad;
489             stepType = NTL_GRADIENT;
490           } else {
491             /* Attempt to use the BFGS direction */
492             PetscCall(MatSolve(tl->M, tao->gradient, tao->stepdirection));
493 
494             /* Check for success (descent direction) */
495             PetscCall(VecDot(tao->stepdirection, tao->gradient, &gdx));
496             if ((gdx <= 0) || PetscIsInfOrNanReal(gdx)) {
497               /* BFGS direction is not descent or direction produced
498                  not a number.  We can assert bfgsUpdates > 1 in this case
499                  Use steepest descent direction (scaled) */
500               PetscCall(MatLMVMReset(tl->M, PETSC_FALSE));
501               PetscCall(MatLMVMUpdate(tl->M, tao->solution, tao->gradient));
502               PetscCall(MatSolve(tl->M, tao->gradient, tao->stepdirection));
503 
504               bfgsUpdates = 1;
505               ++tl->grad;
506               stepType = NTL_GRADIENT;
507             } else {
508               PetscCall(MatLMVMGetUpdateCount(tl->M, &bfgsUpdates));
509               if (1 == bfgsUpdates) {
510                 /* The first BFGS direction is always the scaled gradient */
511                 ++tl->grad;
512                 stepType = NTL_GRADIENT;
513               } else {
514                 ++tl->bfgs;
515                 stepType = NTL_BFGS;
516               }
517             }
518           }
519           break;
520 
521         case NTL_BFGS:
522           /* Can only enter if pc_type == NTL_PC_BFGS
523              Failed to obtain acceptable iterate with BFGS step
524              Attempt to use the scaled gradient direction */
525           PetscCall(MatLMVMReset(tl->M, PETSC_FALSE));
526           PetscCall(MatLMVMUpdate(tl->M, tao->solution, tao->gradient));
527           PetscCall(MatSolve(tl->M, tao->gradient, tao->stepdirection));
528 
529           bfgsUpdates = 1;
530           ++tl->grad;
531           stepType = NTL_GRADIENT;
532           break;
533         }
534         PetscCall(VecScale(tao->stepdirection, -1.0));
535 
536         /* This may be incorrect; linesearch has values for stepmax and stepmin
537            that should be reset. */
538         step = 1.0;
539         PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_reason));
540         PetscCall(TaoAddLineSearchCounts(tao));
541       }
542 
543       if (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER) {
544         /* Failed to find an improving point */
545         f = fold;
546         PetscCall(VecCopy(tl->Xold, tao->solution));
547         PetscCall(VecCopy(tl->Gold, tao->gradient));
548         tao->trust  = 0.0;
549         step        = 0.0;
550         tao->reason = TAO_DIVERGED_LS_FAILURE;
551         break;
552       } else if (stepType == NTL_NEWTON) {
553         if (step < tl->nu1) {
554           /* Very bad step taken; reduce radius */
555           tao->trust = tl->omega1 * PetscMin(norm_d, tao->trust);
556         } else if (step < tl->nu2) {
557           /* Reasonably bad step taken; reduce radius */
558           tao->trust = tl->omega2 * PetscMin(norm_d, tao->trust);
559         } else if (step < tl->nu3) {
560           /* Reasonable step was taken; leave radius alone */
561           if (tl->omega3 < 1.0) {
562             tao->trust = tl->omega3 * PetscMin(norm_d, tao->trust);
563           } else if (tl->omega3 > 1.0) {
564             tao->trust = PetscMax(tl->omega3 * norm_d, tao->trust);
565           }
566         } else if (step < tl->nu4) {
567           /* Full step taken; increase the radius */
568           tao->trust = PetscMax(tl->omega4 * norm_d, tao->trust);
569         } else {
570           /* More than full step taken; increase the radius */
571           tao->trust = PetscMax(tl->omega5 * norm_d, tao->trust);
572         }
573       } else {
574         /* Newton step was not good; reduce the radius */
575         tao->trust = tl->omega1 * PetscMin(norm_d, tao->trust);
576       }
577     } else {
578       /* Trust-region step is accepted */
579       PetscCall(VecCopy(tl->W, tao->solution));
580       f = ftrial;
581       PetscCall(TaoComputeGradient(tao, tao->solution, tao->gradient));
582       ++tl->ntrust;
583     }
584 
585     /* The radius may have been increased; modify if it is too large */
586     tao->trust = PetscMin(tao->trust, tl->max_radius);
587 
588     /* Check for converged */
589     PetscCall(VecNorm(tao->gradient, NORM_2, &gnorm));
590     PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Not-a-Number");
591     needH = 1;
592 
593     PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
594     PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step));
595     PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
596   }
597   PetscFunctionReturn(0);
598 }
599 
600 /* ---------------------------------------------------------- */
601 static PetscErrorCode TaoSetUp_NTL(Tao tao) {
602   TAO_NTL *tl = (TAO_NTL *)tao->data;
603 
604   PetscFunctionBegin;
605   if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient));
606   if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
607   if (!tl->W) PetscCall(VecDuplicate(tao->solution, &tl->W));
608   if (!tl->Xold) PetscCall(VecDuplicate(tao->solution, &tl->Xold));
609   if (!tl->Gold) PetscCall(VecDuplicate(tao->solution, &tl->Gold));
610   tl->bfgs_pre = NULL;
611   tl->M        = NULL;
612   PetscFunctionReturn(0);
613 }
614 
615 /*------------------------------------------------------------*/
616 static PetscErrorCode TaoDestroy_NTL(Tao tao) {
617   TAO_NTL *tl = (TAO_NTL *)tao->data;
618 
619   PetscFunctionBegin;
620   if (tao->setupcalled) {
621     PetscCall(VecDestroy(&tl->W));
622     PetscCall(VecDestroy(&tl->Xold));
623     PetscCall(VecDestroy(&tl->Gold));
624   }
625   PetscCall(KSPDestroy(&tao->ksp));
626   PetscCall(PetscFree(tao->data));
627   PetscFunctionReturn(0);
628 }
629 
630 /*------------------------------------------------------------*/
631 static PetscErrorCode TaoSetFromOptions_NTL(Tao tao, PetscOptionItems *PetscOptionsObject) {
632   TAO_NTL *tl = (TAO_NTL *)tao->data;
633 
634   PetscFunctionBegin;
635   PetscOptionsHeadBegin(PetscOptionsObject, "Newton trust region with line search method for unconstrained optimization");
636   PetscCall(PetscOptionsEList("-tao_ntl_init_type", "radius initialization type", "", NTL_INIT, NTL_INIT_TYPES, NTL_INIT[tl->init_type], &tl->init_type, NULL));
637   PetscCall(PetscOptionsEList("-tao_ntl_update_type", "radius update type", "", NTL_UPDATE, NTL_UPDATE_TYPES, NTL_UPDATE[tl->update_type], &tl->update_type, NULL));
638   PetscCall(PetscOptionsReal("-tao_ntl_eta1", "poor steplength; reduce radius", "", tl->eta1, &tl->eta1, NULL));
639   PetscCall(PetscOptionsReal("-tao_ntl_eta2", "reasonable steplength; leave radius alone", "", tl->eta2, &tl->eta2, NULL));
640   PetscCall(PetscOptionsReal("-tao_ntl_eta3", "good steplength; increase radius", "", tl->eta3, &tl->eta3, NULL));
641   PetscCall(PetscOptionsReal("-tao_ntl_eta4", "excellent steplength; greatly increase radius", "", tl->eta4, &tl->eta4, NULL));
642   PetscCall(PetscOptionsReal("-tao_ntl_alpha1", "", "", tl->alpha1, &tl->alpha1, NULL));
643   PetscCall(PetscOptionsReal("-tao_ntl_alpha2", "", "", tl->alpha2, &tl->alpha2, NULL));
644   PetscCall(PetscOptionsReal("-tao_ntl_alpha3", "", "", tl->alpha3, &tl->alpha3, NULL));
645   PetscCall(PetscOptionsReal("-tao_ntl_alpha4", "", "", tl->alpha4, &tl->alpha4, NULL));
646   PetscCall(PetscOptionsReal("-tao_ntl_alpha5", "", "", tl->alpha5, &tl->alpha5, NULL));
647   PetscCall(PetscOptionsReal("-tao_ntl_nu1", "poor steplength; reduce radius", "", tl->nu1, &tl->nu1, NULL));
648   PetscCall(PetscOptionsReal("-tao_ntl_nu2", "reasonable steplength; leave radius alone", "", tl->nu2, &tl->nu2, NULL));
649   PetscCall(PetscOptionsReal("-tao_ntl_nu3", "good steplength; increase radius", "", tl->nu3, &tl->nu3, NULL));
650   PetscCall(PetscOptionsReal("-tao_ntl_nu4", "excellent steplength; greatly increase radius", "", tl->nu4, &tl->nu4, NULL));
651   PetscCall(PetscOptionsReal("-tao_ntl_omega1", "", "", tl->omega1, &tl->omega1, NULL));
652   PetscCall(PetscOptionsReal("-tao_ntl_omega2", "", "", tl->omega2, &tl->omega2, NULL));
653   PetscCall(PetscOptionsReal("-tao_ntl_omega3", "", "", tl->omega3, &tl->omega3, NULL));
654   PetscCall(PetscOptionsReal("-tao_ntl_omega4", "", "", tl->omega4, &tl->omega4, NULL));
655   PetscCall(PetscOptionsReal("-tao_ntl_omega5", "", "", tl->omega5, &tl->omega5, NULL));
656   PetscCall(PetscOptionsReal("-tao_ntl_mu1_i", "", "", tl->mu1_i, &tl->mu1_i, NULL));
657   PetscCall(PetscOptionsReal("-tao_ntl_mu2_i", "", "", tl->mu2_i, &tl->mu2_i, NULL));
658   PetscCall(PetscOptionsReal("-tao_ntl_gamma1_i", "", "", tl->gamma1_i, &tl->gamma1_i, NULL));
659   PetscCall(PetscOptionsReal("-tao_ntl_gamma2_i", "", "", tl->gamma2_i, &tl->gamma2_i, NULL));
660   PetscCall(PetscOptionsReal("-tao_ntl_gamma3_i", "", "", tl->gamma3_i, &tl->gamma3_i, NULL));
661   PetscCall(PetscOptionsReal("-tao_ntl_gamma4_i", "", "", tl->gamma4_i, &tl->gamma4_i, NULL));
662   PetscCall(PetscOptionsReal("-tao_ntl_theta_i", "", "", tl->theta_i, &tl->theta_i, NULL));
663   PetscCall(PetscOptionsReal("-tao_ntl_mu1", "", "", tl->mu1, &tl->mu1, NULL));
664   PetscCall(PetscOptionsReal("-tao_ntl_mu2", "", "", tl->mu2, &tl->mu2, NULL));
665   PetscCall(PetscOptionsReal("-tao_ntl_gamma1", "", "", tl->gamma1, &tl->gamma1, NULL));
666   PetscCall(PetscOptionsReal("-tao_ntl_gamma2", "", "", tl->gamma2, &tl->gamma2, NULL));
667   PetscCall(PetscOptionsReal("-tao_ntl_gamma3", "", "", tl->gamma3, &tl->gamma3, NULL));
668   PetscCall(PetscOptionsReal("-tao_ntl_gamma4", "", "", tl->gamma4, &tl->gamma4, NULL));
669   PetscCall(PetscOptionsReal("-tao_ntl_theta", "", "", tl->theta, &tl->theta, NULL));
670   PetscCall(PetscOptionsReal("-tao_ntl_min_radius", "lower bound on initial radius", "", tl->min_radius, &tl->min_radius, NULL));
671   PetscCall(PetscOptionsReal("-tao_ntl_max_radius", "upper bound on radius", "", tl->max_radius, &tl->max_radius, NULL));
672   PetscCall(PetscOptionsReal("-tao_ntl_epsilon", "tolerance used when computing actual and predicted reduction", "", tl->epsilon, &tl->epsilon, NULL));
673   PetscOptionsHeadEnd();
674   PetscCall(TaoLineSearchSetFromOptions(tao->linesearch));
675   PetscCall(KSPSetFromOptions(tao->ksp));
676   PetscFunctionReturn(0);
677 }
678 
679 /*------------------------------------------------------------*/
680 static PetscErrorCode TaoView_NTL(Tao tao, PetscViewer viewer) {
681   TAO_NTL  *tl = (TAO_NTL *)tao->data;
682   PetscBool isascii;
683 
684   PetscFunctionBegin;
685   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
686   if (isascii) {
687     PetscCall(PetscViewerASCIIPushTab(viewer));
688     PetscCall(PetscViewerASCIIPrintf(viewer, "Trust-region steps: %" PetscInt_FMT "\n", tl->ntrust));
689     PetscCall(PetscViewerASCIIPrintf(viewer, "Newton search steps: %" PetscInt_FMT "\n", tl->newt));
690     PetscCall(PetscViewerASCIIPrintf(viewer, "BFGS search steps: %" PetscInt_FMT "\n", tl->bfgs));
691     PetscCall(PetscViewerASCIIPrintf(viewer, "Gradient search steps: %" PetscInt_FMT "\n", tl->grad));
692     PetscCall(PetscViewerASCIIPopTab(viewer));
693   }
694   PetscFunctionReturn(0);
695 }
696 
697 /* ---------------------------------------------------------- */
698 /*MC
699   TAONTL - Newton's method with trust region globalization and line search fallback.
700   At each iteration, the Newton trust region method solves the system for d
701   and performs a line search in the d direction:
702 
703             min_d  .5 dT Hk d + gkT d,  s.t.   ||d|| < Delta_k
704 
705   Options Database Keys:
706 + -tao_ntl_init_type - "constant","direction","interpolation"
707 . -tao_ntl_update_type - "reduction","interpolation"
708 . -tao_ntl_min_radius - lower bound on trust region radius
709 . -tao_ntl_max_radius - upper bound on trust region radius
710 . -tao_ntl_epsilon - tolerance for accepting actual / predicted reduction
711 . -tao_ntl_mu1_i - mu1 interpolation init factor
712 . -tao_ntl_mu2_i - mu2 interpolation init factor
713 . -tao_ntl_gamma1_i - gamma1 interpolation init factor
714 . -tao_ntl_gamma2_i - gamma2 interpolation init factor
715 . -tao_ntl_gamma3_i - gamma3 interpolation init factor
716 . -tao_ntl_gamma4_i - gamma4 interpolation init factor
717 . -tao_ntl_theta_i - theta1 interpolation init factor
718 . -tao_ntl_eta1 - eta1 reduction update factor
719 . -tao_ntl_eta2 - eta2 reduction update factor
720 . -tao_ntl_eta3 - eta3 reduction update factor
721 . -tao_ntl_eta4 - eta4 reduction update factor
722 . -tao_ntl_alpha1 - alpha1 reduction update factor
723 . -tao_ntl_alpha2 - alpha2 reduction update factor
724 . -tao_ntl_alpha3 - alpha3 reduction update factor
725 . -tao_ntl_alpha4 - alpha4 reduction update factor
726 . -tao_ntl_alpha4 - alpha4 reduction update factor
727 . -tao_ntl_mu1 - mu1 interpolation update
728 . -tao_ntl_mu2 - mu2 interpolation update
729 . -tao_ntl_gamma1 - gamma1 interpolcation update
730 . -tao_ntl_gamma2 - gamma2 interpolcation update
731 . -tao_ntl_gamma3 - gamma3 interpolcation update
732 . -tao_ntl_gamma4 - gamma4 interpolation update
733 - -tao_ntl_theta - theta1 interpolation update
734 
735   Level: beginner
736 M*/
737 PETSC_EXTERN PetscErrorCode TaoCreate_NTL(Tao tao) {
738   TAO_NTL    *tl;
739   const char *morethuente_type = TAOLINESEARCHMT;
740 
741   PetscFunctionBegin;
742   PetscCall(PetscNewLog(tao, &tl));
743   tao->ops->setup          = TaoSetUp_NTL;
744   tao->ops->solve          = TaoSolve_NTL;
745   tao->ops->view           = TaoView_NTL;
746   tao->ops->setfromoptions = TaoSetFromOptions_NTL;
747   tao->ops->destroy        = TaoDestroy_NTL;
748 
749   /* Override default settings (unless already changed) */
750   if (!tao->max_it_changed) tao->max_it = 50;
751   if (!tao->trust0_changed) tao->trust0 = 100.0;
752 
753   tao->data = (void *)tl;
754 
755   /* Default values for trust-region radius update based on steplength */
756   tl->nu1 = 0.25;
757   tl->nu2 = 0.50;
758   tl->nu3 = 1.00;
759   tl->nu4 = 1.25;
760 
761   tl->omega1 = 0.25;
762   tl->omega2 = 0.50;
763   tl->omega3 = 1.00;
764   tl->omega4 = 2.00;
765   tl->omega5 = 4.00;
766 
767   /* Default values for trust-region radius update based on reduction */
768   tl->eta1 = 1.0e-4;
769   tl->eta2 = 0.25;
770   tl->eta3 = 0.50;
771   tl->eta4 = 0.90;
772 
773   tl->alpha1 = 0.25;
774   tl->alpha2 = 0.50;
775   tl->alpha3 = 1.00;
776   tl->alpha4 = 2.00;
777   tl->alpha5 = 4.00;
778 
779   /* Default values for trust-region radius update based on interpolation */
780   tl->mu1 = 0.10;
781   tl->mu2 = 0.50;
782 
783   tl->gamma1 = 0.25;
784   tl->gamma2 = 0.50;
785   tl->gamma3 = 2.00;
786   tl->gamma4 = 4.00;
787 
788   tl->theta = 0.05;
789 
790   /* Default values for trust region initialization based on interpolation */
791   tl->mu1_i = 0.35;
792   tl->mu2_i = 0.50;
793 
794   tl->gamma1_i = 0.0625;
795   tl->gamma2_i = 0.5;
796   tl->gamma3_i = 2.0;
797   tl->gamma4_i = 5.0;
798 
799   tl->theta_i = 0.25;
800 
801   /* Remaining parameters */
802   tl->min_radius = 1.0e-10;
803   tl->max_radius = 1.0e10;
804   tl->epsilon    = 1.0e-6;
805 
806   tl->init_type   = NTL_INIT_INTERPOLATION;
807   tl->update_type = NTL_UPDATE_REDUCTION;
808 
809   PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
810   PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
811   PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type));
812   PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch, tao));
813   PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix));
814   PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp));
815   PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1));
816   PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix));
817   PetscCall(KSPAppendOptionsPrefix(tao->ksp, "tao_ntl_"));
818   PetscCall(KSPSetType(tao->ksp, KSPSTCG));
819   PetscFunctionReturn(0);
820 }
821