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