xref: /petsc/src/tao/unconstrained/impls/ntr/ntr.c (revision 58d68138c660dfb4e9f5b03334792cd4f2ffd7cc) 
1 #include <../src/tao/unconstrained/impls/ntr/ntrimpl.h>
2 
3 #include <petscksp.h>
4 
5 #define NTR_INIT_CONSTANT      0
6 #define NTR_INIT_DIRECTION     1
7 #define NTR_INIT_INTERPOLATION 2
8 #define NTR_INIT_TYPES         3
9 
10 #define NTR_UPDATE_REDUCTION     0
11 #define NTR_UPDATE_INTERPOLATION 1
12 #define NTR_UPDATE_TYPES         2
13 
14 static const char *NTR_INIT[64] = {"constant", "direction", "interpolation"};
15 
16 static const char *NTR_UPDATE[64] = {"reduction", "interpolation"};
17 
18 /*
19    TaoSolve_NTR - Implements Newton's Method with a trust region approach
20    for solving unconstrained minimization problems.
21 
22    The basic algorithm is taken from MINPACK-2 (dstrn).
23 
24    TaoSolve_NTR computes a local minimizer of a twice differentiable function
25    f by applying a trust region variant of Newton's method.  At each stage
26    of the algorithm, we use the prconditioned conjugate gradient method to
27    determine an approximate minimizer of the quadratic equation
28 
29         q(s) = <s, Hs + g>
30 
31    subject to the trust region constraint
32 
33         || s ||_M <= radius,
34 
35    where radius is the trust region radius and M is a symmetric positive
36    definite matrix (the preconditioner).  Here g is the gradient and H
37    is the Hessian matrix.
38 
39    Note:  TaoSolve_NTR MUST use the iterative solver KSPNASH, KSPSTCG,
40           or KSPGLTR.  Thus, we set KSPNASH, KSPSTCG, or KSPGLTR in this
41           routine regardless of what the user may have previously specified.
42 */
43 static PetscErrorCode TaoSolve_NTR(Tao tao) {
44   TAO_NTR           *tr = (TAO_NTR *)tao->data;
45   KSPType            ksp_type;
46   PetscBool          is_nash, is_stcg, is_gltr, is_bfgs, is_jacobi, is_symmetric, sym_set;
47   KSPConvergedReason ksp_reason;
48   PC                 pc;
49   PetscReal          fmin, ftrial, prered, actred, kappa, sigma, beta;
50   PetscReal          tau, tau_1, tau_2, tau_max, tau_min, max_radius;
51   PetscReal          f, gnorm;
52 
53   PetscReal norm_d;
54   PetscInt  needH;
55 
56   PetscInt i_max = 5;
57   PetscInt j_max = 1;
58   PetscInt i, j, N, n, its;
59 
60   PetscFunctionBegin;
61   if (tao->XL || tao->XU || tao->ops->computebounds) { PetscCall(PetscInfo(tao, "WARNING: Variable bounds have been set but will be ignored by ntr algorithm\n")); }
62 
63   PetscCall(KSPGetType(tao->ksp, &ksp_type));
64   PetscCall(PetscStrcmp(ksp_type, KSPNASH, &is_nash));
65   PetscCall(PetscStrcmp(ksp_type, KSPSTCG, &is_stcg));
66   PetscCall(PetscStrcmp(ksp_type, KSPGLTR, &is_gltr));
67   PetscCheck(is_nash || is_stcg || is_gltr, PetscObjectComm((PetscObject)tao), PETSC_ERR_SUP, "TAO_NTR requires nash, stcg, or gltr for the KSP");
68 
69   /* Initialize the radius and modify if it is too large or small */
70   tao->trust = tao->trust0;
71   tao->trust = PetscMax(tao->trust, tr->min_radius);
72   tao->trust = PetscMin(tao->trust, tr->max_radius);
73 
74   /* Allocate the vectors needed for the BFGS approximation */
75   PetscCall(KSPGetPC(tao->ksp, &pc));
76   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &is_bfgs));
77   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &is_jacobi));
78   if (is_bfgs) {
79     tr->bfgs_pre = pc;
80     PetscCall(PCLMVMGetMatLMVM(tr->bfgs_pre, &tr->M));
81     PetscCall(VecGetLocalSize(tao->solution, &n));
82     PetscCall(VecGetSize(tao->solution, &N));
83     PetscCall(MatSetSizes(tr->M, n, n, N, N));
84     PetscCall(MatLMVMAllocate(tr->M, tao->solution, tao->gradient));
85     PetscCall(MatIsSymmetricKnown(tr->M, &sym_set, &is_symmetric));
86     PetscCheck(sym_set && is_symmetric, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix in the LMVM preconditioner must be symmetric.");
87   } else if (is_jacobi) PetscCall(PCJacobiSetUseAbs(pc, PETSC_TRUE));
88 
89   /* Check convergence criteria */
90   PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
91   PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm));
92   PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Inf or NaN");
93   needH = 1;
94 
95   tao->reason = TAO_CONTINUE_ITERATING;
96   PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
97   PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, 1.0));
98   PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
99   if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0);
100 
101   /*  Initialize trust-region radius */
102   switch (tr->init_type) {
103   case NTR_INIT_CONSTANT:
104     /*  Use the initial radius specified */
105     break;
106 
107   case NTR_INIT_INTERPOLATION:
108     /*  Use the initial radius specified */
109     max_radius = 0.0;
110 
111     for (j = 0; j < j_max; ++j) {
112       fmin  = f;
113       sigma = 0.0;
114 
115       if (needH) {
116         PetscCall(TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre));
117         needH = 0;
118       }
119 
120       for (i = 0; i < i_max; ++i) {
121         PetscCall(VecCopy(tao->solution, tr->W));
122         PetscCall(VecAXPY(tr->W, -tao->trust / gnorm, tao->gradient));
123         PetscCall(TaoComputeObjective(tao, tr->W, &ftrial));
124 
125         if (PetscIsInfOrNanReal(ftrial)) {
126           tau = tr->gamma1_i;
127         } else {
128           if (ftrial < fmin) {
129             fmin  = ftrial;
130             sigma = -tao->trust / gnorm;
131           }
132 
133           PetscCall(MatMult(tao->hessian, tao->gradient, tao->stepdirection));
134           PetscCall(VecDot(tao->gradient, tao->stepdirection, &prered));
135 
136           prered = tao->trust * (gnorm - 0.5 * tao->trust * prered / (gnorm * gnorm));
137           actred = f - ftrial;
138           if ((PetscAbsScalar(actred) <= tr->epsilon) && (PetscAbsScalar(prered) <= tr->epsilon)) {
139             kappa = 1.0;
140           } else {
141             kappa = actred / prered;
142           }
143 
144           tau_1   = tr->theta_i * gnorm * tao->trust / (tr->theta_i * gnorm * tao->trust + (1.0 - tr->theta_i) * prered - actred);
145           tau_2   = tr->theta_i * gnorm * tao->trust / (tr->theta_i * gnorm * tao->trust - (1.0 + tr->theta_i) * prered + actred);
146           tau_min = PetscMin(tau_1, tau_2);
147           tau_max = PetscMax(tau_1, tau_2);
148 
149           if (PetscAbsScalar(kappa - (PetscReal)1.0) <= tr->mu1_i) {
150             /*  Great agreement */
151             max_radius = PetscMax(max_radius, tao->trust);
152 
153             if (tau_max < 1.0) {
154               tau = tr->gamma3_i;
155             } else if (tau_max > tr->gamma4_i) {
156               tau = tr->gamma4_i;
157             } else {
158               tau = tau_max;
159             }
160           } else if (PetscAbsScalar(kappa - (PetscReal)1.0) <= tr->mu2_i) {
161             /*  Good agreement */
162             max_radius = PetscMax(max_radius, tao->trust);
163 
164             if (tau_max < tr->gamma2_i) {
165               tau = tr->gamma2_i;
166             } else if (tau_max > tr->gamma3_i) {
167               tau = tr->gamma3_i;
168             } else {
169               tau = tau_max;
170             }
171           } else {
172             /*  Not good agreement */
173             if (tau_min > 1.0) {
174               tau = tr->gamma2_i;
175             } else if (tau_max < tr->gamma1_i) {
176               tau = tr->gamma1_i;
177             } else if ((tau_min < tr->gamma1_i) && (tau_max >= 1.0)) {
178               tau = tr->gamma1_i;
179             } else if ((tau_1 >= tr->gamma1_i) && (tau_1 < 1.0) && ((tau_2 < tr->gamma1_i) || (tau_2 >= 1.0))) {
180               tau = tau_1;
181             } else if ((tau_2 >= tr->gamma1_i) && (tau_2 < 1.0) && ((tau_1 < tr->gamma1_i) || (tau_2 >= 1.0))) {
182               tau = tau_2;
183             } else {
184               tau = tau_max;
185             }
186           }
187         }
188         tao->trust = tau * tao->trust;
189       }
190 
191       if (fmin < f) {
192         f = fmin;
193         PetscCall(VecAXPY(tao->solution, sigma, tao->gradient));
194         PetscCall(TaoComputeGradient(tao, tao->solution, tao->gradient));
195 
196         PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm));
197         PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Inf or NaN");
198         needH = 1;
199 
200         PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
201         PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, 1.0));
202         PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
203         if (tao->reason != TAO_CONTINUE_ITERATING) { PetscFunctionReturn(0); }
204       }
205     }
206     tao->trust = PetscMax(tao->trust, max_radius);
207 
208     /*  Modify the radius if it is too large or small */
209     tao->trust = PetscMax(tao->trust, tr->min_radius);
210     tao->trust = PetscMin(tao->trust, tr->max_radius);
211     break;
212 
213   default:
214     /*  Norm of the first direction will initialize radius */
215     tao->trust = 0.0;
216     break;
217   }
218 
219   /* Have not converged; continue with Newton method */
220   while (tao->reason == TAO_CONTINUE_ITERATING) {
221     /* Call general purpose update function */
222     PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
223     ++tao->niter;
224     tao->ksp_its = 0;
225     /* Compute the Hessian */
226     if (needH) {
227       PetscCall(TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre));
228       needH = 0;
229     }
230 
231     if (tr->bfgs_pre) {
232       /* Update the limited memory preconditioner */
233       PetscCall(MatLMVMUpdate(tr->M, tao->solution, tao->gradient));
234     }
235 
236     while (tao->reason == TAO_CONTINUE_ITERATING) {
237       PetscCall(KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre));
238 
239       /* Solve the trust region subproblem */
240       PetscCall(KSPCGSetRadius(tao->ksp, tao->trust));
241       PetscCall(KSPSolve(tao->ksp, tao->gradient, tao->stepdirection));
242       PetscCall(KSPGetIterationNumber(tao->ksp, &its));
243       tao->ksp_its += its;
244       tao->ksp_tot_its += its;
245       PetscCall(KSPCGGetNormD(tao->ksp, &norm_d));
246 
247       if (0.0 == tao->trust) {
248         /* Radius was uninitialized; use the norm of the direction */
249         if (norm_d > 0.0) {
250           tao->trust = norm_d;
251 
252           /* Modify the radius if it is too large or small */
253           tao->trust = PetscMax(tao->trust, tr->min_radius);
254           tao->trust = PetscMin(tao->trust, tr->max_radius);
255         } else {
256           /* The direction was bad; set radius to default value and re-solve
257              the trust-region subproblem to get a direction */
258           tao->trust = tao->trust0;
259 
260           /* Modify the radius if it is too large or small */
261           tao->trust = PetscMax(tao->trust, tr->min_radius);
262           tao->trust = PetscMin(tao->trust, tr->max_radius);
263 
264           PetscCall(KSPCGSetRadius(tao->ksp, tao->trust));
265           PetscCall(KSPSolve(tao->ksp, tao->gradient, tao->stepdirection));
266           PetscCall(KSPGetIterationNumber(tao->ksp, &its));
267           tao->ksp_its += its;
268           tao->ksp_tot_its += its;
269           PetscCall(KSPCGGetNormD(tao->ksp, &norm_d));
270 
271           PetscCheck(norm_d != 0.0, PetscObjectComm((PetscObject)tao), PETSC_ERR_PLIB, "Initial direction zero");
272         }
273       }
274       PetscCall(VecScale(tao->stepdirection, -1.0));
275       PetscCall(KSPGetConvergedReason(tao->ksp, &ksp_reason));
276       if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) && (tr->bfgs_pre)) {
277         /* Preconditioner is numerically indefinite; reset the
278            approximate if using BFGS preconditioning. */
279         PetscCall(MatLMVMReset(tr->M, PETSC_FALSE));
280         PetscCall(MatLMVMUpdate(tr->M, tao->solution, tao->gradient));
281       }
282 
283       if (NTR_UPDATE_REDUCTION == tr->update_type) {
284         /* Get predicted reduction */
285         PetscCall(KSPCGGetObjFcn(tao->ksp, &prered));
286         if (prered >= 0.0) {
287           /* The predicted reduction has the wrong sign.  This cannot
288              happen in infinite precision arithmetic.  Step should
289              be rejected! */
290           tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d);
291         } else {
292           /* Compute trial step and function value */
293           PetscCall(VecCopy(tao->solution, tr->W));
294           PetscCall(VecAXPY(tr->W, 1.0, tao->stepdirection));
295           PetscCall(TaoComputeObjective(tao, tr->W, &ftrial));
296 
297           if (PetscIsInfOrNanReal(ftrial)) {
298             tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d);
299           } else {
300             /* Compute and actual reduction */
301             actred = f - ftrial;
302             prered = -prered;
303             if ((PetscAbsScalar(actred) <= tr->epsilon) && (PetscAbsScalar(prered) <= tr->epsilon)) {
304               kappa = 1.0;
305             } else {
306               kappa = actred / prered;
307             }
308 
309             /* Accept or reject the step and update radius */
310             if (kappa < tr->eta1) {
311               /* Reject the step */
312               tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d);
313             } else {
314               /* Accept the step */
315               if (kappa < tr->eta2) {
316                 /* Marginal bad step */
317                 tao->trust = tr->alpha2 * PetscMin(tao->trust, norm_d);
318               } else if (kappa < tr->eta3) {
319                 /* Reasonable step */
320                 tao->trust = tr->alpha3 * tao->trust;
321               } else if (kappa < tr->eta4) {
322                 /* Good step */
323                 tao->trust = PetscMax(tr->alpha4 * norm_d, tao->trust);
324               } else {
325                 /* Very good step */
326                 tao->trust = PetscMax(tr->alpha5 * norm_d, tao->trust);
327               }
328               break;
329             }
330           }
331         }
332       } else {
333         /* Get predicted reduction */
334         PetscCall(KSPCGGetObjFcn(tao->ksp, &prered));
335         if (prered >= 0.0) {
336           /* The predicted reduction has the wrong sign.  This cannot
337              happen in infinite precision arithmetic.  Step should
338              be rejected! */
339           tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
340         } else {
341           PetscCall(VecCopy(tao->solution, tr->W));
342           PetscCall(VecAXPY(tr->W, 1.0, tao->stepdirection));
343           PetscCall(TaoComputeObjective(tao, tr->W, &ftrial));
344           if (PetscIsInfOrNanReal(ftrial)) {
345             tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
346           } else {
347             PetscCall(VecDot(tao->gradient, tao->stepdirection, &beta));
348             actred = f - ftrial;
349             prered = -prered;
350             if ((PetscAbsScalar(actred) <= tr->epsilon) && (PetscAbsScalar(prered) <= tr->epsilon)) {
351               kappa = 1.0;
352             } else {
353               kappa = actred / prered;
354             }
355 
356             tau_1   = tr->theta * beta / (tr->theta * beta - (1.0 - tr->theta) * prered + actred);
357             tau_2   = tr->theta * beta / (tr->theta * beta + (1.0 + tr->theta) * prered - actred);
358             tau_min = PetscMin(tau_1, tau_2);
359             tau_max = PetscMax(tau_1, tau_2);
360 
361             if (kappa >= 1.0 - tr->mu1) {
362               /* Great agreement; accept step and update radius */
363               if (tau_max < 1.0) {
364                 tao->trust = PetscMax(tao->trust, tr->gamma3 * norm_d);
365               } else if (tau_max > tr->gamma4) {
366                 tao->trust = PetscMax(tao->trust, tr->gamma4 * norm_d);
367               } else {
368                 tao->trust = PetscMax(tao->trust, tau_max * norm_d);
369               }
370               break;
371             } else if (kappa >= 1.0 - tr->mu2) {
372               /* Good agreement */
373 
374               if (tau_max < tr->gamma2) {
375                 tao->trust = tr->gamma2 * PetscMin(tao->trust, norm_d);
376               } else if (tau_max > tr->gamma3) {
377                 tao->trust = PetscMax(tao->trust, tr->gamma3 * norm_d);
378               } else if (tau_max < 1.0) {
379                 tao->trust = tau_max * PetscMin(tao->trust, norm_d);
380               } else {
381                 tao->trust = PetscMax(tao->trust, tau_max * norm_d);
382               }
383               break;
384             } else {
385               /* Not good agreement */
386               if (tau_min > 1.0) {
387                 tao->trust = tr->gamma2 * PetscMin(tao->trust, norm_d);
388               } else if (tau_max < tr->gamma1) {
389                 tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
390               } else if ((tau_min < tr->gamma1) && (tau_max >= 1.0)) {
391                 tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
392               } else if ((tau_1 >= tr->gamma1) && (tau_1 < 1.0) && ((tau_2 < tr->gamma1) || (tau_2 >= 1.0))) {
393                 tao->trust = tau_1 * PetscMin(tao->trust, norm_d);
394               } else if ((tau_2 >= tr->gamma1) && (tau_2 < 1.0) && ((tau_1 < tr->gamma1) || (tau_2 >= 1.0))) {
395                 tao->trust = tau_2 * PetscMin(tao->trust, norm_d);
396               } else {
397                 tao->trust = tau_max * PetscMin(tao->trust, norm_d);
398               }
399             }
400           }
401         }
402       }
403 
404       /* The step computed was not good and the radius was decreased.
405          Monitor the radius to terminate. */
406       PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
407       PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, tao->trust));
408       PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
409     }
410 
411     /* The radius may have been increased; modify if it is too large */
412     tao->trust = PetscMin(tao->trust, tr->max_radius);
413 
414     if (tao->reason == TAO_CONTINUE_ITERATING) {
415       PetscCall(VecCopy(tr->W, tao->solution));
416       f = ftrial;
417       PetscCall(TaoComputeGradient(tao, tao->solution, tao->gradient));
418       PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm));
419       PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Inf or NaN");
420       needH = 1;
421       PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
422       PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, tao->trust));
423       PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
424     }
425   }
426   PetscFunctionReturn(0);
427 }
428 
429 /*------------------------------------------------------------*/
430 static PetscErrorCode TaoSetUp_NTR(Tao tao) {
431   TAO_NTR *tr = (TAO_NTR *)tao->data;
432 
433   PetscFunctionBegin;
434   if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient));
435   if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
436   if (!tr->W) PetscCall(VecDuplicate(tao->solution, &tr->W));
437 
438   tr->bfgs_pre = NULL;
439   tr->M        = NULL;
440   PetscFunctionReturn(0);
441 }
442 
443 /*------------------------------------------------------------*/
444 static PetscErrorCode TaoDestroy_NTR(Tao tao) {
445   TAO_NTR *tr = (TAO_NTR *)tao->data;
446 
447   PetscFunctionBegin;
448   if (tao->setupcalled) { PetscCall(VecDestroy(&tr->W)); }
449   PetscCall(KSPDestroy(&tao->ksp));
450   PetscCall(PetscFree(tao->data));
451   PetscFunctionReturn(0);
452 }
453 
454 /*------------------------------------------------------------*/
455 static PetscErrorCode TaoSetFromOptions_NTR(Tao tao, PetscOptionItems *PetscOptionsObject) {
456   TAO_NTR *tr = (TAO_NTR *)tao->data;
457 
458   PetscFunctionBegin;
459   PetscOptionsHeadBegin(PetscOptionsObject, "Newton trust region method for unconstrained optimization");
460   PetscCall(PetscOptionsEList("-tao_ntr_init_type", "tao->trust initialization type", "", NTR_INIT, NTR_INIT_TYPES, NTR_INIT[tr->init_type], &tr->init_type, NULL));
461   PetscCall(PetscOptionsEList("-tao_ntr_update_type", "radius update type", "", NTR_UPDATE, NTR_UPDATE_TYPES, NTR_UPDATE[tr->update_type], &tr->update_type, NULL));
462   PetscCall(PetscOptionsReal("-tao_ntr_eta1", "step is unsuccessful if actual reduction < eta1 * predicted reduction", "", tr->eta1, &tr->eta1, NULL));
463   PetscCall(PetscOptionsReal("-tao_ntr_eta2", "", "", tr->eta2, &tr->eta2, NULL));
464   PetscCall(PetscOptionsReal("-tao_ntr_eta3", "", "", tr->eta3, &tr->eta3, NULL));
465   PetscCall(PetscOptionsReal("-tao_ntr_eta4", "", "", tr->eta4, &tr->eta4, NULL));
466   PetscCall(PetscOptionsReal("-tao_ntr_alpha1", "", "", tr->alpha1, &tr->alpha1, NULL));
467   PetscCall(PetscOptionsReal("-tao_ntr_alpha2", "", "", tr->alpha2, &tr->alpha2, NULL));
468   PetscCall(PetscOptionsReal("-tao_ntr_alpha3", "", "", tr->alpha3, &tr->alpha3, NULL));
469   PetscCall(PetscOptionsReal("-tao_ntr_alpha4", "", "", tr->alpha4, &tr->alpha4, NULL));
470   PetscCall(PetscOptionsReal("-tao_ntr_alpha5", "", "", tr->alpha5, &tr->alpha5, NULL));
471   PetscCall(PetscOptionsReal("-tao_ntr_mu1", "", "", tr->mu1, &tr->mu1, NULL));
472   PetscCall(PetscOptionsReal("-tao_ntr_mu2", "", "", tr->mu2, &tr->mu2, NULL));
473   PetscCall(PetscOptionsReal("-tao_ntr_gamma1", "", "", tr->gamma1, &tr->gamma1, NULL));
474   PetscCall(PetscOptionsReal("-tao_ntr_gamma2", "", "", tr->gamma2, &tr->gamma2, NULL));
475   PetscCall(PetscOptionsReal("-tao_ntr_gamma3", "", "", tr->gamma3, &tr->gamma3, NULL));
476   PetscCall(PetscOptionsReal("-tao_ntr_gamma4", "", "", tr->gamma4, &tr->gamma4, NULL));
477   PetscCall(PetscOptionsReal("-tao_ntr_theta", "", "", tr->theta, &tr->theta, NULL));
478   PetscCall(PetscOptionsReal("-tao_ntr_mu1_i", "", "", tr->mu1_i, &tr->mu1_i, NULL));
479   PetscCall(PetscOptionsReal("-tao_ntr_mu2_i", "", "", tr->mu2_i, &tr->mu2_i, NULL));
480   PetscCall(PetscOptionsReal("-tao_ntr_gamma1_i", "", "", tr->gamma1_i, &tr->gamma1_i, NULL));
481   PetscCall(PetscOptionsReal("-tao_ntr_gamma2_i", "", "", tr->gamma2_i, &tr->gamma2_i, NULL));
482   PetscCall(PetscOptionsReal("-tao_ntr_gamma3_i", "", "", tr->gamma3_i, &tr->gamma3_i, NULL));
483   PetscCall(PetscOptionsReal("-tao_ntr_gamma4_i", "", "", tr->gamma4_i, &tr->gamma4_i, NULL));
484   PetscCall(PetscOptionsReal("-tao_ntr_theta_i", "", "", tr->theta_i, &tr->theta_i, NULL));
485   PetscCall(PetscOptionsReal("-tao_ntr_min_radius", "lower bound on initial trust-region radius", "", tr->min_radius, &tr->min_radius, NULL));
486   PetscCall(PetscOptionsReal("-tao_ntr_max_radius", "upper bound on trust-region radius", "", tr->max_radius, &tr->max_radius, NULL));
487   PetscCall(PetscOptionsReal("-tao_ntr_epsilon", "tolerance used when computing actual and predicted reduction", "", tr->epsilon, &tr->epsilon, NULL));
488   PetscOptionsHeadEnd();
489   PetscCall(KSPSetFromOptions(tao->ksp));
490   PetscFunctionReturn(0);
491 }
492 
493 /*------------------------------------------------------------*/
494 /*MC
495   TAONTR - Newton's method with trust region for unconstrained minimization.
496   At each iteration, the Newton trust region method solves the system.
497   NTR expects a KSP solver with a trust region radius.
498             min_d  .5 dT Hk d + gkT d,  s.t.   ||d|| < Delta_k
499 
500   Options Database Keys:
501 + -tao_ntr_init_type - "constant","direction","interpolation"
502 . -tao_ntr_update_type - "reduction","interpolation"
503 . -tao_ntr_min_radius - lower bound on trust region radius
504 . -tao_ntr_max_radius - upper bound on trust region radius
505 . -tao_ntr_epsilon - tolerance for accepting actual / predicted reduction
506 . -tao_ntr_mu1_i - mu1 interpolation init factor
507 . -tao_ntr_mu2_i - mu2 interpolation init factor
508 . -tao_ntr_gamma1_i - gamma1 interpolation init factor
509 . -tao_ntr_gamma2_i - gamma2 interpolation init factor
510 . -tao_ntr_gamma3_i - gamma3 interpolation init factor
511 . -tao_ntr_gamma4_i - gamma4 interpolation init factor
512 . -tao_ntr_theta_i - theta1 interpolation init factor
513 . -tao_ntr_eta1 - eta1 reduction update factor
514 . -tao_ntr_eta2 - eta2 reduction update factor
515 . -tao_ntr_eta3 - eta3 reduction update factor
516 . -tao_ntr_eta4 - eta4 reduction update factor
517 . -tao_ntr_alpha1 - alpha1 reduction update factor
518 . -tao_ntr_alpha2 - alpha2 reduction update factor
519 . -tao_ntr_alpha3 - alpha3 reduction update factor
520 . -tao_ntr_alpha4 - alpha4 reduction update factor
521 . -tao_ntr_alpha4 - alpha4 reduction update factor
522 . -tao_ntr_mu1 - mu1 interpolation update
523 . -tao_ntr_mu2 - mu2 interpolation update
524 . -tao_ntr_gamma1 - gamma1 interpolcation update
525 . -tao_ntr_gamma2 - gamma2 interpolcation update
526 . -tao_ntr_gamma3 - gamma3 interpolcation update
527 . -tao_ntr_gamma4 - gamma4 interpolation update
528 - -tao_ntr_theta - theta interpolation update
529 
530   Level: beginner
531 M*/
532 
533 PETSC_EXTERN PetscErrorCode TaoCreate_NTR(Tao tao) {
534   TAO_NTR *tr;
535 
536   PetscFunctionBegin;
537 
538   PetscCall(PetscNewLog(tao, &tr));
539 
540   tao->ops->setup          = TaoSetUp_NTR;
541   tao->ops->solve          = TaoSolve_NTR;
542   tao->ops->setfromoptions = TaoSetFromOptions_NTR;
543   tao->ops->destroy        = TaoDestroy_NTR;
544 
545   /* Override default settings (unless already changed) */
546   if (!tao->max_it_changed) tao->max_it = 50;
547   if (!tao->trust0_changed) tao->trust0 = 100.0;
548   tao->data = (void *)tr;
549 
550   /*  Standard trust region update parameters */
551   tr->eta1 = 1.0e-4;
552   tr->eta2 = 0.25;
553   tr->eta3 = 0.50;
554   tr->eta4 = 0.90;
555 
556   tr->alpha1 = 0.25;
557   tr->alpha2 = 0.50;
558   tr->alpha3 = 1.00;
559   tr->alpha4 = 2.00;
560   tr->alpha5 = 4.00;
561 
562   /*  Interpolation trust region update parameters */
563   tr->mu1 = 0.10;
564   tr->mu2 = 0.50;
565 
566   tr->gamma1 = 0.25;
567   tr->gamma2 = 0.50;
568   tr->gamma3 = 2.00;
569   tr->gamma4 = 4.00;
570 
571   tr->theta = 0.05;
572 
573   /*  Interpolation parameters for initialization */
574   tr->mu1_i = 0.35;
575   tr->mu2_i = 0.50;
576 
577   tr->gamma1_i = 0.0625;
578   tr->gamma2_i = 0.50;
579   tr->gamma3_i = 2.00;
580   tr->gamma4_i = 5.00;
581 
582   tr->theta_i = 0.25;
583 
584   tr->min_radius = 1.0e-10;
585   tr->max_radius = 1.0e10;
586   tr->epsilon    = 1.0e-6;
587 
588   tr->init_type   = NTR_INIT_INTERPOLATION;
589   tr->update_type = NTR_UPDATE_REDUCTION;
590 
591   /* Set linear solver to default for trust region */
592   PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp));
593   PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1));
594   PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix));
595   PetscCall(KSPAppendOptionsPrefix(tao->ksp, "tao_ntr_"));
596   PetscCall(KSPSetType(tao->ksp, KSPSTCG));
597   PetscFunctionReturn(0);
598 }
599